The first microwave and submillimetre closure study using particle models of oriented ice hydrometeors to simulate polarimetric measurements of ice clouds

The first closure study involving passive microwave and submillimetre measurements of ice clouds with the consideration of oriented particles is presented, using a unique combination of polarised observations from the ISMAR spectral-like radiometer, two radars with frequencies of 35 and 95 GHz, and a variety of in situ instruments. Of particular interest to this study are the large V–H polarised brightness temperature differences measured from ISMAR above a thick frontal ice cloud. Previous studies combining radar and passive submillimetre measurements have not considered polarisation differences. Moreover, they have assumed particle habits a priori. We aim to test whether the large V–H measurements can be simulated successfully by using an atmospheric model consistent with in situ microphysics.

An atmospheric model is constructed using information from the in situ measurements, such as the ice water content, the particle size distribution, and the mass and shape of particles, as well as background information obtained from dropsonde profiles. Columnar and dendritic aggregate particle models are generated specifically for this case, and their scattering properties are calculated using the independent monomer approximation under the assumption of horizontal orientation. The scattering properties are used to perform polarised radiative transfer simulations using ARTS to test whether we can successfully simulate the measured large V–H differences. Radar measurements are used to extrapolate the 1-D microphysical profile to derive a time series of particle size distributions which are used to simulate ISMAR brightness temperatures. These simulations are compared to the observations.

It is found that particle models that are consistent with in situ microphysics observations are capable of reproducing the brightness temperature depression and polarisation signature measured from ISMAR at the dual-polarised channel of 243 GHz. However, it was required that a proportion of the particles were changed in order to increase the V–H polarised brightness temperature differences. Thus, we incorporated millimetre-sized dendritic crystals, as these particles were observed in the probe imagery. At the second dual-polarised channel of 664 GHz, the brightness temperature depressions were generally simulated at the correct locations; however, the simulated V–H was too large. This work shows that multi-frequency polarisation information could be used to infer realistic particle shapes, orientations, and representations of the split between single crystals and aggregates within the cloud.

1 Introduction

Passive radiometry allows for measurements of the column mass of atmospheric ice, since millimetre and submillimetre waves are sensitive to scattering by ice. Efforts have been made to improve airborne and spaceborne retrievals of ice water path (IWP) by measuring submillimetre brightness temperatures (e.g. Evans et al., 2005, 2012; Brath et al., 2018; Fox et al., 2017; Kangas et al., 2014). The Ice Cloud Imager (ICI) will be the first operational instrument to cover submillimetre wavelengths, with frequencies ranging from 183 to 664 GHz (Eriksson et al., 2020). The instrument has been specifically designed for measuring cloud ice from space and is due for launch on board the MetOp-SG satellite “B” in 2025. It is expected that the combination of frequencies available on ICI will allow for more accurate estimations of IWP and mean mass dimension, as there is a dependence between these properties and the submillimetre brightness temperature depression (Evans et al., 1998, 2002; Buehler et al., 2007). However, the retrieval of IWP using submillimetre-wave measurements will depend on the microphysical and particle size distribution (PSD) assumptions (Baran et al., 2018; Fox et al., 2019). Previous studies have found that by comparing the brightness temperatures at simultaneous orthogonal horizontal and vertical polarisations, which in this paper we denote by H and V, it is possible to gain some information about the size, shape, and orientation of ice particles within the cloud (e.g. Evans and Stephens, 1995; Miao et al., 2003; Xie and Miao, 2011; Defer et al., 2014; Ding et al., 2017; Gong and Wu, 2017; Zhang and Gasiewski, 2018). For that reason, two channels on ICI have the capability of measuring at both horizontal and vertical polarisations, namely 243 and 664 GHz. Indeed, to prepare for polarised brightness temperature measurements from ICI, recent studies have worked towards representing hydrometeor orientation in a simplified manner in data assimilation and retrieval applications (Barlakas et al., 2021; Kaur et al., 2022).

Covering a frequency range of 118 to 874 GHz, the International Sub-Millimetre Airborne Radiometer (ISMAR) has been developed by the Met Office and ESA as an airborne demonstrator instrument for ICI. The instrument flies on the Facility for Airborne Atmospheric Measurements (FAAM) BAe-146 research aircraft and is useful for testing ice-scattering models that could be used within retrieval algorithms for ICI. There have already been successful applications of the ISMAR in retrieving IWP and surface emissivity properties by Brath et al. (2018) and Prigent et al. (2017), respectively. More recently, Fox et al. (2019) presented a microwave and submillimetre closure study which utilised the Microwave Airborne Radiometer Scanning System (MARSS) (McGrath and Hewison, 2001) and ISMAR observations between 183 and 664 GHz obtained from above a few cases of mid-latitude cirrus. These measurements were concurrently simulated using the Atmospheric Radiative Transfer Simulator (ARTS; Eriksson et al., 2011; Buehler et al., 2018) with a number of ice crystal models from the corresponding single-scattering database (Eriksson et al., 2018). The masses of the particles were constrained by the in situ bulk ice water content (IWC) measurements. Even this closure study, where the details of the in situ shapes, sizes, and bulk IWC were accounted for, showed that no single ice crystal model from the Eriksson et al. (2018) database could fully replicate the observations at all the frequencies considered simultaneously, and in one case, none of the assumed models replicated the measurements at all the considered frequencies. The study of Fox et al. (2019) indicates the difficulties that might be encountered in utilising ICI measurements to retrieve IWP. Fox (2020) used output from a high-resolution NWP model in radiative transfer simulations and showed that the simulations are sensitive to assumed particle shape, particularly at 243 GHz. However, they pointed out that in a case with greater ice mass at higher altitudes, stronger sensitivity at higher frequencies would be expected. The above studies did not consider polarisation differences.

In this study, a comprehensive microwave closure experiment has been performed, using data collected during flight B984 of the North Atlantic Waveguide and Downstream Impact EXperiment (NAWDEX) campaign. This flight was also studied by Ewald et al. (2021) and Pfreundschuh et al. (2022). The authors of the latter study focussed on retrieving IWC and vertical distributions of ice hydrometeors. However, they assumed particle habits a priori, while our study is unique as we choose ice particle models that are consistent with the in situ cloud measurements. Moreover, we focus on forward-modelling brightness temperatures and polarisation differences, while polarisation was not considered by the aforementioned authors. Measurements of deep frontal cloud were obtained in a region off the west coast of Scotland on 14 October 2016. Three aircraft collected coincident measurements from above the cloud during this case, at an altitude of approximately 9.5 km. These are listed in Sect. 2.1 and described in more detail in Schäfler et al. (2018). Independent datasets from a variety of in situ and remote-sensing instruments, such as ISMAR, are utilised. The overall goal is to determine whether we can combine the in situ measurements to construct an atmospheric model that can successfully replicate the large brightness temperature depressions and V–H differences measured from ISMAR. Thus, we perform polarised radiative transfer simulations using ARTS. The differences in brightness temperatures between those simulated using ARTS and those measured by the airborne radiometer ISMAR are analysed. To quantify the importance of utilising polarisation observations at microwave and submillimetre frequencies, we focus mainly on ISMAR polarisation measurements at 243 GHz, along with comparing a number of simulations and measurements at 664 GHz.

The choice of how to represent cloud ice and snow in radiative transfer models has developed greatly over the past 2 decades. As particle size increases with respect to the wavelength, the particle shape and structure play a significant role in different interference patterns that are found within the crystal (e.g. McCusker et al., 2019; Kleanthous et al., 2022). Waves scattered by different parts of the ice particle interfere in the far field, which affects the far-field scattering properties. Thus, the far-field scattering is sensitive to the structure of the ice particles on scales comparable to the wavelength of the radiation, and as measurements move into the millimetre and submillimetre part of the spectrum, the scattering properties become increasingly sensitive to the detailed structure of the ice particles. Thus, the community has shifted from employing drastic simplifications of cloud ice and snow (i.e. approximating particles by spheres or spheroids of equivalent size) to using a more realistic representation of ice particles. For example, the current version of RTTOV-SCATT (v13.0) uses a large plate aggregate (Geer et al., 2021). For the simulations performed here, we tune the choice of ice particle ensemble by using in situ measurements. We use cloud particle imagery from the CIP cloud imaging probes to choose particle habits for the simulations. Mass–size relationships are derived specifically for this case, and particles are generated to match the relationships. This process is described in more detail in the following sections. We combine radar reflectivity measurements at 35 and 95 GHz in conjunction with ISMAR polarisation differences (V–H) to test if the in situ-derived ice aggregate models can consistently replicate these observations across the microwave spectrum. Furthermore, we consider the added value of utilising radar reflectivity observations in tandem with passive microwave dual-polarisation observations. The instrumentation used and measurements obtained from both above and within the cloud are described in more detail below.

2 Details of the case study

2.1 Above-cloud measurements

Figure 1a shows the region of the flight path of the aircraft, as well as the surface rain rate estimated from an operational Met Office C-band radar. The FAAM aircraft carried the ISMAR radiometer, which measured brightness temperatures at an off-nadir observation angle of between 51 and 52°, matching the planned configuration of ICI. A bias correction was applied to the vertically polarised 664 GHz data, as described by Fox (2020). Coincident data were also obtained from two different radars. The 35 GHz radar was on board the German High Altitude and LOng Range Research Aircraft (HALO) and is part of the HALO Microwave Package (HAMP; Mech et al., 2014). The 95 GHz RAdar SysTem Airborne (RASTA) Doppler cloud radar, part of the RALI synergistic radar–lidar platform, flew on board the French Service des Avions Français Instrumentés pour la Recherche en Environnement (SAFIRE) Falcon 20 aircraft. More information on the RASTA radar can be found in Delanoë et al. (2013), while comparisons between HAMP and RASTA can be found in Ewald et al. (2019). Both radars operated at an observation angle near 0°, i.e. nadir.

Figure 1b shows the altitude and latitude of the FAAM aircraft. The ISMAR brightness temperatures were measured during the above-cloud run, between approximately 10:00 and 10:20 UTC. In situ cloud measurements were obtained following the above-cloud run, between 10:37 and 11:00 UTC. More detail on the in situ measurements is found in Sect. 2.2.

Figure 1 (a) Rain rate in mm h⁻¹ estimated from the Met Office C-band radar (5.6 GHz) at Druim a'Starraig in Scotland at 10:15 UTC. The red star shows the location of the radar. A black and yellow dashed line has been overlaid between latitudes of approximately 56.5–59.5° N at a longitude of 6.5° W to highlight the region of analysis where the three aircraft flew during the above-cloud near-coincident run. (b) Altitude and latitude of the FAAM aircraft during the time of interest. The aircraft performed an above-cloud run, during which the ISMAR brightness temperatures were measured. Following that, in situ cloud measurements were obtained during a profile descent.

Figure 2 shows the radar reflectivities measured at (a) 35 GHz and (b) 95 GHz, and crosses are shown above the cloud tops at latitudes where dropsondes were released. It is obvious that the 35 GHz radar with its high peak pulse power of 30 kW is more sensitive, detecting a signal from particles in the upper region of the cloud that are not picked up by the higher-frequency radar (which has a lower power of 1.8 kW). Each panel has been divided into four regions, labelled a–d above the panels. These represent different cloud regions as follows: a – small-scale convective cloud with precipitation, along with overlying mid- and upper-level cloud layers; b – broken up mid- and high-level cloud visible at 35 GHz; c – deep frontal cloud with intermittent precipitation; d – thinning frontal cloud with intermittent precipitation. Two regions of particularly high reflectivity are seen in region c, where the 35 and 95 GHz radars measure values of about 20 and 10 dBZ, respectively. These are located at altitudes between about 2 and 4 km and latitudes between approximately 58.4 and 57.9°, as well as 57.7 and 57.3°. A thin melting layer is observed at ≈ 1.5 km, with a weak bright band evident in the 35 GHz radar reflectivity.

Figure 2 Measured reflectivities from (a) HAMP (35 GHz) and (b) RASTA (95 GHz). Each radar image has been divided into four parts and labelled a–d, representing different cloud regions as follows: a – small-scale convective cloud with precipitation, along with overlying mid- and upper-level cloud layers; b – broken up mid- and high-level cloud visible at 35 GHz; c – deep frontal cloud with intermittent precipitation; d – thinning frontal cloud with intermittent precipitation. Crosses are shown above the cloud tops at latitudes where dropsondes were released. Panel (c) shows the brightness temperatures at H and V polarisations measured at different latitudes using the ISMAR radiometer at 243 GHz, and panel (d) shows the V–H brightness temperature difference.

Figure 2c shows the 243 GHz brightness temperatures at H and V polarisations measured from ISMAR during this campaign, and Fig. 2d shows the V–H polarimetric differences. The shaded regions represent estimates of the measurement uncertainties, described in more detail in Appendix A.

Figure 2c shows that large brightness temperature depressions were measured as the aircraft flew over the deep frontal cloud region (region c of the radar data), with a decrease of approximately 30 K at a latitude of 57.3°. Figure 2d shows that there were also regions where a large V–H polarimetric signal was measured, reaching almost 10 K. The V–H signal is correlated with the brightness temperature depression, indicating that it is microphysical rather than being caused by the surface. Moreover, the high reflectivities in Fig. 2a and b are correlated with the large brightness temperature depressions and V–H measurements and thus could be caused by large oriented ice particles. Comparison of the radar and ISMAR data in Fig. 2 shows that there are regions of cloud without an obvious melting layer but a large V–H difference, for example at latitudes close to 57.6–57.5°. This is evidence that the large polarimetric signal is caused by the ice cloud particles above the melting layer. Thus, although the polarisation difference is generally increased as a result of the melting layer (e.g. Gong and Wu, 2017), we do not attempt to represent scattering by the melting layer in our simulations.

We note that there is a positional offset between the ISMAR and radar measurements (corresponding to approximately 0.1° latitude at the ground) due to the forward-viewing setup of ISMAR. This is discussed further later in the paper.

2.2 In situ measurements

As mentioned above, we utilise in situ measurements in order to perform the radiative transfer simulations. In particular, we model a multi-layer atmosphere using in situ PSDs along with particle models generated specifically for this case. This is described in more detail in Sect. 3. The in situ measurements were obtained following the above-cloud run, when the FAAM aircraft turned around and performed a profile descent along the same track (see Fig. 1b). Thus, all the observations were not obtained coincidently in time, resulting in a limitation of the experiment. This is discussed further in Sect. 5, where we develop a technique to overcome the limitation by using radar reflectivities to derive PSDs to use at other times.

The aircraft carried the CIP-15 and CIP-100 optical array probes (OAPs), measuring particles of 15 to 960 µm, as well as 100 to 6400 µm, respectively. The probes provide measurements of the PSDs, along with 2-D imagery of particles, which is useful to decide which particle habit to use for the simulations. A review of the different OAPs and their characteristics is given by McFarquhar et al. (2017). Note that for any given particle the maximum dimension measured from the 2-D image will be less than or equal to the true maximum dimension. However, for a large enough sample size the maximum size measured in the PSD will be close to the maximum dimension of the sampled particles, even for randomly oriented particles (because at least some of the particles will have their true maximum dimension oriented close to the imaging plane). Also, since the images are measured in the horizontal plane, any tendency for the particles to be horizontally oriented will reduce the difference between the maximum size estimated from the 2-D image and the true maximum dimension. To obtain a PSD covering a wide range of particle sizes, the PSDs from the two CIP instruments were composited using the method described in Cotton et al. (2013). The deep cone Nevzorov probe provided data on the liquid and total water (ice plus liquid) contents. The probe provides accurate IWC measurements for ice particles with $D_{max}<$ 4 mm (Korolev et al., 2013a), but there are larger measurement uncertainties for sizes > 4 mm due to fragmentation and shattering. The measured PSDs and IWC were averaged within each layer.

3 Simulation setup and construction of a model atmosphere

ARTS is particularly useful for this work because of its capability to handle polarised radiative transfer. Brightness temperatures are calculated using the RT4 polarised radiative transfer model within ARTS, described in Evans and Stephens (1995), and a viewing angle of 50° is used. Construction of an atmospheric model is required in order to perform the radiative transfer simulations. Prior to descending through the cloud, the aircraft released a series of dropsondes to obtain the background atmospheric state for input into ARTS. The latitudes at which these were released are shown by black crosses in Fig. 2. As well as providing water vapour mixing ratio profiles of the atmosphere, the dropsonde profiles provide surface properties such as temperature and wind speed that feed into the TESSEM ocean surface emissivity model (Prigent et al., 2017) used within ARTS. We note that the aircraft is over land at some points, roughly between 58.4–57.9, 57.5–57.4, and 56.7–56.6°. We discuss the impact of surface emissivity on brightness temperatures in Sect. 4.1. The Rosenkranz (1998) gas absorption model is used in ARTS.

As a supplement to the sounding, it is necessary to input information on the ice cloud that was present during the study. Due to the depth of the cloud, it would not be possible to obtain an adequate representation of the atmospheric conditions using a single homogeneous layer. Hence, we represent vertical variation of the PSDs and particle shapes by modelling the atmosphere using seven different layers with depths of approximately 1 km each. These layers are located between altitudes of 2 and 9 km. RT4 assumes that particles are azimuthally random; i.e. there is some preferential polar alignment, but the particles are randomly oriented in the azimuth. Furthermore, the model assumes a plane-parallel atmosphere. Thus, in all the simulations performed here, we use a 1-D plane-parallel atmosphere within ARTS, with azimuthally randomly oriented particles (more information on the particles is found in Sect. 3.1). Barlakas and Eriksson (2020) investigated the impact of 3-D effects on millimetre and submillimetre brightness temperatures measured by satellite radiometers. They found that the difference between approximating the 3-D scene by a 1-D plane-parallel approximation was dominated by the heterogeneity of the cloud field within the beam (which is small for airborne instruments such as ISMAR) and that horizontal photon transport between different parts of the scene was small (brightness temperature differences < 1 K typically). Although their study did not consider polarimetric effects, it seems likely that the uncertainties in V–H will be of comparably small magnitude, and we argue that 3-D effects can be reasonably neglected in our study, in the context of other uncertainties (in particular, the lack of colocation between the in situ sampling and the radiometer measurements). However, we acknowledge that more research on the influence of 3-D radiative transfer effects on polarised brightness temperatures is required in order to fully quantify this. Precipitation was intermittently present below 2 km, with a thin melting layer showing as a bright band in the radar reflectivity in Fig. 2a and rain below it. As mentioned, we do not represent melting particles here and assume a Marshall–Palmer distribution of rain beneath the ice cloud base. The distribution used here corresponds to a rain rate of 1–2 mm h⁻¹, estimated from the radar data shown in Fig. 1a.

In addition to including a database of scattering calculations for realistic particle habits performed using the discrete dipole approximation (DDA; Eriksson et al., 2018), ARTS also accepts externally generated scattering calculations. This gives us the opportunity to generate particles specifically for this study, along with providing the flexibility to use an alternative scattering method for the calculations. These points are described in more detail in the following subsections.

3.1 Particle generation

To decide on which particle habits to use, imagery from the CIP cloud imaging probes is used. Examples of imagery from different altitudes are shown in Fig. 3a. Panel (a) shows a slide from each of the seven cloud layers, imaged from the CIP-15 probe. Also shown on the right side of panel (a) are examples of dendritic particles selected from the CIP-100 imagery in the bottom layer of cloud.

Figure 3 (a) CIP-15 images from each of the seven cloud layers profiled by the aircraft. The height of each frame is approximately 0.96 mm. The layer altitudes are given to the right of the particle imagery. We approximate the top three layers (labelled L1–L3 in this figure) as columnar aggregates and the bottom four layers (L4–L7) as dendritic aggregates. Examples of the generated model particles are shown on the left side of the panel. Examples of dendritic particles imaged in L7 by the CIP-100 probe are also shown on the right side of the panel. (b) The final mass–size relationships used to model the particles for this study, using SI units. Also plotted are relationships derived by Brown and Francis (1995) and Cotton et al. (2013).

Mixtures of particle habits were present throughout the cloud, but visual inspection of the imagery led us to approximate the atmospheric model using two different particle types. We use columnar aggregates higher up in the cloud, between 9 and 6 km, and dendritic aggregates lower in the cloud, between 6 and 2 km (i.e. layers labelled L4–L7 in Fig. 3a). We note that dendritic monomers of 1–2 mm in size were also present in L7, as can be seen from the CIP-100 imagery. An ensemble of ice particles over a range of different sizes were generated for each cloud layer. The particles follow a realistic mass–size relationship that is consistent with the observed cloud properties. The process of choosing the mass–size relationships is explained further in Sect. 3.2. A range of columnar and dendritic aggregates were constructed using the particle aggregation model of Westbrook et al. (2004). Within the model, the monomer shape and size were specified at the outset, and realistic aggregates were generated via the mechanism of differential sedimentation. Those with masses that match the derived mass–size relationships to within 20 % were kept, storing a range of sizes to represent the full PSD. Examples of the generated particles are shown on the left side of Fig. 3a. The monomers are oriented randomly upon generation, as are the resulting aggregates. We assume the particles are horizontally oriented but are randomly oriented in the horizontal plane (i.e. azimuthally random orientation) to match the assumptions made in RT4. To do this, the particles generated by the aggregation model are reoriented based on their maximum moment of inertia, such that the maximum distribution of mass is in the horizontal plane, following the same method used by Tyynelä et al. (2011). In this article we refer to such particles as “horizontally oriented”. We do not assume a particular canting angle. The aspect ratio of the generated aggregates is calculated as the ratio of the maximum particle size along the z axis to the maximum size in the x–y plane.

The particle size bins corresponding to the distribution suggest that the cloud contained particles up to approximately $D_{max}=$ 5.75 mm. There are large uncertainties in the number concentrations of ice particles smaller than 100 µm, due to shattering (Korolev et al., 2013b). Furthermore, Buehler et al. (2007) show the sensitivity of various submillimetre channels to particles of different size. They show that for a fixed IWP, brightness temperature differences are not sensitive to particles smaller than 100 µm at the lower frequencies considered here, while at higher frequencies close to 664 GHz, there is slightly more (but still very little) sensitivity to these small particles. Thus, we do not include particles smaller than 100 µm here. As described above, a population of particles was generated using the aggregation model. For sizes larger than $D_{max}=$ 100 µm, the population was subsampled to obtain a range of particles that are representative of the sizes measured in each layer (i.e. the size bins corresponding to the layer-averaged composited PSDs), while the smallest particles were ignored by setting the concentrations of any measured particles below 100 µm to 0. The final number of aggregate realisations used for each layer ranges from 46 to 62. Although we did not perform a detailed sensitivity test for the number of realisations used, we note that our sampling strategy is comparable to Eriksson et al. (2018), who used approximately 35–45 particles to span across the range of sizes relevant to ice hydrometeors.

3.2 Derivation of mass–size relationships

The relationship between particle mass and maximum dimension, D_max, is usually assumed to have the form $m=a{D}_{max}^b$, where m is the mass of the particle, and the parameters a and b are constants which depend on particle habit and atmospheric conditions such as temperature. For example, Cotton et al. (2013) proposed $m=\mathrm{0.0257}{D}_{max}^{\mathrm{2.12}}$ as a good fit to several ice clouds. As discussed by Mason et al. (2018) and references therein, the prefactor, a, of the mass–size relationship scales the ice particle's density, and the exponent, b, is related to the particle shape or growth mechanism. Mitchell et al. (1990) presented mass–size relationships with values of b ranging from 1.7 to 2.6 for different habits. As discussed by Westbrook et al. (2004), the value of b for unrimed aggregate snowflakes is usually around 2. Note that mass-dimension parameterisations consider D_max to be the maximum dimension of the particle in any direction. This is slightly different compared to OAP measurements, in that the instruments measure particles in two perpendicular directions, with D_max commonly calculated by fitting the smallest circle to fully enclose the image.

For the simulations performed here, we aim to do as much as possible to try and match the atmospheric state at the time of the in situ observations. Therefore, rather than employing relationships that are commonly used in the literature, a realistic mass–size relationship is derived for each layer. We use two measurements to constrain our choice of mass–size relationship for each layer, namely the ice water content measured using the Nevzorov probe, as well as the radar reflectivity, Z. The IWC is measured at the same time and location as the PSD, thus providing a direct constraint on possible combinations of a and b in the m–D_max relationships. However, as mentioned in Sect. 2.2, there are large measurement uncertainties when particles are > 4 mm. Z is measured at a different time and location, thus providing an additional but weaker constraint. Therefore, our first step to find suitable combinations of a and b for this case is to simulate the IWC with various possible m–D_max relationships and compare to the IWC measured using the Nevzorov probe. We start from the Cotton et al. (2013) relationship of $m=\mathrm{0.0257}{D}_{max}^{\mathrm{2.12}}$ as a baseline. First we fix a= 0.0257 kg m^−b and vary b until good agreement is found with the measured bulk IWC data (differences within 1 %), and then we fix b=2.12 and vary a to match the measured IWC. This leaves us with two sets of potential relationships to choose from. We then use Z to refine our choice. We generate aggregates (as described in Sect. 3.1) to match each of the two sets of relationships and use the independent monomer approximation (IMA) to calculate the radar cross section, σ_r, of the generated particles at 35 and 95 GHz. These calculations are used to simulate the above-cloud equivalent radar reflectivity, Z_e, in order to test the suitability of the generated particles for this case. The equation for Z_e is given by Atlas et al. (1995)

$$\begin{array}{}\text{(1)}& Z_{\mathrm{e}}={\mathrm{10}}^{\mathrm{18}}C\int {\mathit{\sigma }}_{\mathrm{r}}\left(D_{max}\right)n\left(D_{max}\right)\mathrm{d}{D}_{max},\end{array}$$

where $C={\mathit{\lambda }}^{\mathrm{4}}/\left({\mathit{\pi }}^{\mathrm{5}}|\left({\mathit{ϵ}}_{\mathrm{liquid}}-\mathrm{1}\right)/\left({\mathit{ϵ}}_{\mathrm{liquid}}+\mathrm{2}\right){|}^{\mathrm{2}}\right)$ is a frequency-dependent constant, and n(D_max) represents the in situ distribution of particles. Multiplication by 10¹⁸ converts the units of Z_e from m³ to conventional radar meteorology units of mm⁶ m⁻³. More information on the method can be found in Baran et al. (2014). Two-way attenuation by ice was estimated for this case and determined to be negligible, with values ≪ 1 dBZ_e at 35 and 95 GHz. Thus, we do not include it in any of the calculations presented here.

After comparing simulated to measured Z_e, the final m–D_max relationship for each layer is then chosen by selecting the one whose simulated Z_e lies closest to the mean of the measured Z_e distribution. In most layers, the exponent b has a value of 2.12, except in layer 2 where better results were found when a was fixed to 0.0257 kg m^−b. The mass–size relationships are plotted in Fig. 3b, along with commonly used versions of Brown and Francis (1995) and Cotton et al. (2013). The relationships derived for the lower layers of cloud are very consistent with the aforementioned relationships, while particles in the top layers of cloud have lower masses that would not be represented correctly by those of Brown and Francis (1995) and Cotton et al. (2013). As expected, b≈ 2 at the bottom of the cloud where large aggregate snowflakes are present.

Figure 4 shows 2-D histograms of the measured reflectivities, with the colour bar representing the number of observations in each bin. Circles are overlaid at the central altitude of each of the seven ice cloud layers, showing the simulated reflectivity using the final modelled particles and the layer-averaged PSDs. Black lines are plotted at the left side of the reflectivities, showing the estimated sensitivity of each of the radars. Below this noise level, the signal is not detectable by the radar. The minimum detectable signal is calculated in dBZ as 10log ₁₀(r²)+c, where the range, r, is equivalent to the aircraft altitude minus height, and c is a constant. Simulations of the 35 GHz reflectivities in Fig. 4a are in good agreement with the measurements. The simulated value in the bottom layer is slightly larger than the main bulk of the measured distribution, though it still falls within the observed values. In the case of the 95 GHz simulations in Fig. 4b, the low sensitivity of the radar means there is a very clear line in the measurements below which no signal is picked up, and simulations in the top three layers are below the minimum detectable signal. However, we show in Sect. 4.1 that the ISMAR brightness temperatures and polarisation signal at 243 GHz are not sensitive to these top layers of cloud anyway. In the following section, we use these particles to simulate the H and V polarisation measurements from ISMAR.

Figure 4 Two-dimensional histograms of observed reflectivities from the (a) HAMP 35 GHz and (b) RASTA 95 GHz radars, with the colour bar representing the number of observations in each bin. Simulated reflectivities for each layer are plotted with circles on top of the measurements. The circles show the reflectivities calculated using the layer-averaged PSDs along with the horizontally aligned particles generated to follow the mass–size relationships derived from measurements. Black lines at the left edge of the reflectivities have been plotted to show the estimated noise level, below which no signal is detected.

3.3 Particle scattering

In this study we calculate the scattering properties of ice aggregates using the independent monomer approximation (IMA), outlined in McCusker et al. (2020). Like DDA (e.g. Draine and Flatau, 1994, and references therein), the method involves discretising a particle into volume elements which are treated as dipoles. Whereas the DDA method considers interactions between all dipoles, IMA only considers interaction between dipoles that are in the same monomer. This simplifies calculations and allows considerable reductions to time and memory requirements, particularly for aggregates with a large number of monomers. Defining the size parameter as $x=k{D}_{max}/\mathrm{2}$, where $k=\mathrm{2}\mathit{\pi }/\mathit{\lambda }$ is the wavenumber corresponding to an incident wave of wavelength λ, McCusker et al. (2020) showed that scattering calculations using the method closely agree with DDA for x less than ∼ 5, with biases of less than 10 % in the scattering cross sections for all monomer habits considered in that study (i.e. plates, columns, and dendrites). A size parameter of x= 5 corresponds to quite small particle sizes of $D_{max}=$ 2 mm at 243 GHz and 0.7 mm at 664 GHz. However, aggregates comprising dendritic monomers have significantly lower biases, remaining within 10 % at considerably larger values of x up to 18 (Tyynelä et al., 2023), corresponding to $D_{max}=$ 7 mm at 243 GHz and 2.6 mm at 664 GHz. The aforementioned studies also show that IMA can successfully reproduce radar multi-wavelength and multi-polarisation parameters.

To verify the accuracy of IMA for this application, tests were performed by simulating brightness temperatures at 243 GHz using a monodispersive distribution of particles within each layer, with the mean particle size in each layer chosen to match PSD measurements. Results using three different scattering approximations were compared, specifically IMA, DDA, and RGA (i.e. the Rayleigh–Gans approximation; see Bohren and Huffman, 1983). The brightness temperatures simulated using DDA were within 0.05 K of the IMA results, while RGA failed to produce an equivalent depression, differing from the DDA results by up to 2.2 K. With RGA the brightness temperature depressions would be much smaller than the more accurate IMA technique, thereby resulting in erroneous simulations.

In the atmospheric model used here, we approximate the top three cloud layers with columnar aggregates. As mentioned, McCusker et al. (2020) showed that IMA generally performs better for fluffier dendritic aggregates compared to aggregates with more compact monomers such as columns and plates. The largest columnar particle generated for layer 3 (we show later in Fig. 5 that simulations are highly sensitive to this layer at 664 GHz but less so at 243 GHz) is 2.5 mm, corresponding to x≈ 17. This is larger than the maximum size parameter of 10 tested for columnar aggregates by McCusker et al. (2020). Thus, the columnar shape combined with the large size parameters may lead to scattering errors at 664 GHz. However, we performed tests (not shown here) which revealed that the strongest contribution to scattering at 664 GHz in all layers comes from small particles less than 2 mm. Comparisons with DDA calculations revealed that IMA is accurate to within 10 % for these particles. Thus, we are confident that the IMA method is sufficiently accurate to be applied to this study at both 243 and 664 GHz.

Yurkin and Hoekstra (2007) suggest 10 dipoles per wavelength as a rule of thumb for discretisation of particles in the DDA, and given that IMA is a simplification of DDA, we apply this same rule of thumb to our analysis. All particles considered at 243 GHz had > 10 dipoles per wavelength. As mentioned, scattering in all layers is dominated by sub-2 mm particles at 664 GHz, and all particles < 2 mm fulfil the criteria of 10 dipoles per wavelength at 664 GHz. While some of our largest aggregates were more coarsely sampled at this high frequency, our analysis is not sensitive to the precise representation of the largest particles. Scattering is simulated using a grid with 5° spacing for the incident and scattering polar angles and the scattering azimuth angle. Orientation-averaged scattering is calculated using 36 azimuthal orientations (using a regular azimuth angle grid with 10° spacing). This is sufficient for the particles used here, with tests showing that the mean error in the first phase matrix element using 36 orientations is within 0.4 % of the results using 360 orientations (1° azimuth angle grid). The scattering properties are integrated over the ensemble of particles weighted by the particle size distribution.

4 Simulation of the ISMAR polarised brightness temperatures at 243 GHz

In this section, we focus mainly on 243 GHz. We can estimate the size parameters (x) of interest at this frequency from the particle measurements. In the top three cloud layers, which we are modelling with columnar aggregates, particle sizes up to approximately 2.5 mm were measured, corresponding to x=6.4. Larger particles were measured in the lower cloud layers (which we model with dendritic aggregates), corresponding to x= 13.8. Considering the discussion above, we estimate the bias in IMA scattering cross sections to be within about 10 % in both cases. Although our main focus is on 243 GHz, we elaborate further on the results at 664 GHz in Sect. 6. As described in the previous section, a range of particles have been generated for each of the seven model layers to follow the in situ measurements as closely as possible. Here, simulations are performed using the full range of generated particle models along with the measured PSDs.

4.1 Multi-layer, polydispersive distribution

The results for the polydispersive case are shown in Fig. 5. The simulations have been performed using the PSDs for each of the seven layers, shown by the crosses in Fig. 6a, with each cross representing one of the generated particles. The single-scattering properties of each of the particles are calculated using IMA and incorporated into ARTS along with the PSDs. One layer of cloud has been added at a time, until the full model cloud is included. This provides insight into which parts of the cloud profile weight the observed brightness temperature depressions most. In Fig. 5, the points along the abscissa represent the gradual increase in cloud layers used in the simulation, starting with the clear-sky case. Then a Marshall–Palmer distribution of rain is inserted between the ground and the cloud base at 2 km. The values used are $N_{\mathrm{0}}=\mathrm{8}\times {\mathrm{10}}^{\mathrm{6}}$ m⁻⁴ and $\mathit{\lambda }=\mathrm{4}\times {\mathrm{10}}^{\mathrm{3}}$ m⁻¹, which corresponds to a precipitation rate of 1.12 mm h⁻¹ with the Marshall Palmer distribution. This rain distribution is included in all further simulations presented here.

Figure 5 (a, c) Brightness temperature and (b, d) V–H brightness temperature difference at 243 GHz (a, b) and 664 GHz (c, d). The first point along the abscissa is the clear-sky case, and the second point represents the case where the model atmosphere includes a rain distribution but no cloud. The third until final points show results when one layer of cloud is added to the model atmosphere (with rain) at a time, starting from the top layer of cloud only (third point) and ending with all seven cloud layers (final point).

Figure 6 (a) The lines with crosses show the PSD number concentrations, and the straight lines show exponential fits to the measurements. The exponential fits are introduced in Sect. 5. Note that the minimum value of D_max is 100 µm. (b) Two-dimensional histogram of the ISMAR V–H data, with the colour bar representing the number of observations in each bin. Markers have been overlaid to show the simulation results. The red star shows the V–H brightness temperature difference calculated using the measured PSDs (i.e. the results in Fig. 5a and b when all layers of the model atmosphere are included), and the yellow diamonds show the values simulated when the aggregates in the lowest region of cloud are replaced with horizontally aligned single dendrites. The different diamonds show results obtained when the amount of cloud comprising dendrites is gradually increased, as described in the text.

The third point along the x axis displays the result when the top layer of ice cloud is included, along with the rain distribution. Then the second layer of ice is added, and so on until the full seven layers of cloud ice along with a distribution of rain below the cloud base are included. The plot is done in this way to mimic the increasing depressions measured by ISMAR as the aircraft flew over the cloud, as in Fig. 2c. The individual H and V brightness temperatures at 243 GHz are displayed in Fig. 5a, with V–H shown in Fig. 5b. V and H are depressed to values around 240 K. Note that the equivalent calculations were also performed using RGA (not shown here). It was found that the RGA scattering method underestimates brightness temperature depressions that can be simulated using IMA, by up to 9.6 K.

At 243 GHz, there is a polarisation signal of about 2 K in the clear-sky case, which comes from the surface emissivity model. Including rain, represented by liquid spheres, diminishes the polarisation signal, meaning the modelled signal is due to the cloud only. The top three layers of cloud (9–6 km) have little effect on brightness temperature depressions or polarisation. As additional model layers are included in the simulation, the brightness temperatures become more depressed at both H and V. The depression increases due to the signal from layers 4–7 (6–2 km), particularly layer 7. However, layers 4–6 (6–3 km) have the greatest effect on the polarisation signal (0.3–1.5 K), with a small decrease with inclusion of layer 7. It is interesting to note that when layer 7 is included, the brightness temperature depression increases further but the polarisation signal decreases slightly to below 3 K. This is perhaps a result of including too many large, lower-density particles in layer 7, which have a weak polarisation dependence. McCusker et al. (2020) show that the IMA can successfully reproduce polarimetric parameters such as the differential reflectivity (up to 200 GHz). Since IMA only includes interactions within individual monomers, this implies that the monomer shapes and distribution of monomer orientations within the aggregate determine the polarisation properties rather than the shape of the “envelope” around the aggregate. Thus, a decrease in simulated V–H may be due to the aspect ratio or orientation of the monomer crystals within the aggregates in layer 7. It is possible that the dendritic monomers should be oriented before aggregation or pivoted on attachment, both of which result in flatter and more dense particles (Schrom et al., 2023). Alternatively, it may be that the aspect ratio of the monomers needs to be smaller (i.e. thinner dendrites).

The same experiment of adding layers gradually was also performed at 664 GHz (Fig. 5c and d). At that frequency there is no polarisation signal from the surface. The brightness temperature depression is affected by the mid-high cloud layers 1–4, particularly 3 and 4 (7–5 km), while the depression does not increase further with inclusion of lower cloud in layers 5–7. The polarisation signal is also mainly sensitive to layers 3 and 4 (1.9, 3.4 K), but shows some sensitivity to layers 2 and 5 (0.6–0.7 K) and a very small amount of sensitivity to layers 1 and 6 (0.1–0.2 K). The bottom cloud layer does not increase the polarisation signal.

Because of the measurement strategy adopted for this case study, we only have a single microphysical profile, and therefore we simulate only one value of brightness temperature and V–H. However, we have a time series of these variables measured from ISMAR. Since the cloud is heterogeneous, it is not reasonable to compare the simulation against specific samples from ISMAR. Instead, we look at whether the simulated (V, V–H) lies within the distribution of the measured samples. Comparisons of the simulated V–H at 243 GHz (using the complete model atmosphere) with the ISMAR measurements are shown in Fig. 6b. The ISMAR measurements are plotted as a 2-D histogram, and the simulated V–H brightness temperature difference is shown by the red star. Although the simulation falls within the ISMAR measurements, the V–H value of approximately 3 K is relatively low for such a large brightness temperature depression, considering values of up to 8 K were measured. It is interesting to consider why the simulated polarimetric signal is not as great as that measured by ISMAR and explore ways in which it could be increased to match the observations. We noted above that the addition of layer 7 in our model atmosphere causes a decrease in simulated V–H. As mentioned above, the small polarisation signal could be caused by the inclusion of too many large, lower-density particles in the bottom layer of cloud, which do not have a strong polarisation signature, or the aspect ratio or orientation of the monomer crystals within the aggregates may not be realistic enough. The particle habit used in that layer of our model may not adequately represent the real cloud. A potential reason for this could be that there was a change in microphysics between the times of the ISMAR measurements which were made at 10:00–10:20 UTC and the in situ cloud measurements taken between 10:37 and 11:00 UTC. Gong and Wu (2017) showed that the V–H differences can be increased by changing the particle habit or aspect ratio. We investigate these points in the following subsection. It is worth noting that a further way of increasing the polarimetric difference would be to consider habit mixture models, such as by Miao et al. (2003). Using aggregates of differing monomer shapes in our simulations may increase the V–H differences, but this is beyond the scope of the present study.

4.2 Changing particle habit

A potential reason for the small polarimetric signal in the original simulations is that aggregates are not responsible for the V–H brightness temperature difference observed from ISMAR. Up until this point, monomers have not been included in our model atmosphere, and there is evidence from the imagery in Fig. 3a that such particles were indeed present between 2 and 3 km at the time of interest. Thus, in this subsection we incorporate horizontally aligned dendrites into the lowest layer of cloud. It is worth pointing out that since we are using single crystals, the scattering calculations for these additional dendrites are done using DDA rather than IMA, as the IMA method is only applicable to aggregates.

A monodispersive distribution of horizontally aligned single dendrites is used to replace the aggregates in the lowest portion of the cloud. The dendrites have a size of approximately $D_{max}=$ 1 mm, similar to what was imaged by the CIP-100 probe (Fig. 3a). Calculating the aspect ratio as the ratio between the length of the particle in the z direction and the maximum width of the particle in the x–y plane (i.e. values closer to 1 are more spherical while smaller values indicate a thinner dendrite), the single dendrites have an aspect ratio of approximately 0.1. The number concentration is chosen in such a way that the measured IWC over the 1 km-deep layer is maintained. Different heights of the cloud layer are replaced by dendrites, starting with the bottom 100 m and increasing the height by 100 m at a time until finally the lowest 500 m of cloud is replaced with dendrites. In other words, we replace 10 %–50 % of the bottom layer of cloud. These values correspond to dendrites comprising approximately 3 % to 15 % of the total IWP.

Figure 6b shows the V–H brightness temperature differences plotted along with the values measured from ISMAR. The red star shows the original result before the inclusion of dendrites. The smallest diamond closest to the star is the result when 10 % of the bottom layer (i.e. 100 m of cloud) is replaced with dendrites, and the largest, lowest diamond shows the result obtained when 50 % of the bottom layer (i.e. the lowest 500 m of cloud) is replaced. Even when no dendrites are included, the largest brightness temperature depression is almost captured. However, adding dendrites increases the polarisation difference, with each extra 100 m increasing V–H by approximately 1–1.5 K, while V remains almost constant. A V–H value of 8.5 K is obtained when the lowest 400 m of the model cloud is replaced with single dendrites, while V–H reaches 10.1 K when the lowest 500 m is replaced. These values are very consistent with ISMAR measurements, showing that agreement with observations is possible by including a small IWP of oriented dendritic ice crystals in the simulation.

5 Simulating a time series of ISMAR measurements at 243 GHz using synthetic PSDs that accurately reproduce Z at 35 GHz

One of the limitations of the available measurements is the fact that the in situ profile is not coincident with the brightness temperature measurements. Moreover, using layer-averaged PSDs for each of the seven cloud layers results in simulating only one value of brightness temperature. Clearly this does not represent the heterogeneity of the cloud. Here we try to represent a time series of brightness temperatures, analogous to what was measured by ISMAR. This is done using exponential PSDs $N\left(D_{max}\right)=N_{\mathrm{0}}\exp \left(-\mathit{\lambda }{D}_{max}\right)$, which were obtained from the measured PSDs in each layer by plotting D_max against ln (N) and fitting straight lines to the distribution (i.e. the straight lines in Fig. 6a). For consistency with the original setup, the parameterised PSDs have been truncated using the lower and upper limits of the measured PSDs. We then adjust the exponential PSDs that were fit to the layer-averaged PSDs. Since Gong and Wu (2017) showed that V–H can be increased by increasing the mean size of the PSD, we choose to adjust the PSDs by fixing the N₀ values while changing λ. One could alternatively fix λ and change N₀ (i.e. change the number density of particles); however, this would not increase the mean particle size.

Data from the 35 GHz radar have been used for this work. As mentioned, the radar was on board the HALO aircraft, while the ISMAR was on board the FAAM aircraft. Since the two aircraft were travelling at different speeds (with HALO travelling approximately 60 m s⁻¹ faster than FAAM), it is not possible to directly compare time series data, since different parts of the cloud were measured at different times. Because of this, we bin the data from both platforms by latitude, assuming that the cloud scene at each point in space does not change substantially in the interval between the radar measurement and the ISMAR measurement. The 35 GHz radar data were sorted into 200 latitude bins, and an average reflectivity, $\stackrel{\mathrm{‾}}{Z}$, for each bin was calculated for each of the seven layers. For each value of $\stackrel{\mathrm{‾}}{Z}$, the exponential PSDs are adjusted by changing λ such that a simulation of the radar reflectivity using our generated particles matches $\stackrel{\mathrm{‾}}{Z}$ at that latitude. Thus, a synthetic time series of in situ PSDs is generated for each layer, which is then loaded into ARTS to simulate the brightness temperature at different latitudes, or equivalently different times.

One aspect we need to consider here is whether the model particle ensemble generated for this study is sufficient for simulations performed with different PSDs. By plotting the integrand of the reflectivity (N(D)σ(D)) using different PSDs, we may determine whether larger particles would contribute significantly to the total scattering of the ensemble. If N(D)σ(D) has a bell shape with a clear peak and tails, the generated model particles are sufficient. Otherwise, if N(D)σ(D) is truncated, the PSD contains particles which are too big to be represented by the model particle ensemble. The straight line that was fit to the measured PSDs in the bottom layer (layer 7) has a λ value of 1301 m⁻¹, while the λ values calculated for the time series range from 353 to 21 800 m⁻¹. Here we find that values of λ less than about 1700 m⁻¹ result in a distribution that is too broad to be fully represented using the particles that were generated for this case; i.e. the results would be affected by larger particles in the distribution than what we have generated here. This means we cannot trust our retrievals of λ in these cases, so we set a threshold value of λ> 1700 m⁻¹. Alternatively, one could generate larger particles to represent the full distribution. However, we choose not to do that here since larger particles were not measured by the in situ instruments.

As mentioned previously, there is imperfect collocation of the ISMAR and radar measurements since the radar views at nadir and ISMAR is forward-looking, viewing at ≈ 51°. This means a correction must be applied to account for the positional offset between ISMAR and the radar in the time series simulations. This offset varies with height, but since the simulations at 243 GHz are most sensitive to the bottom layers of cloud, we will consider cloud between 2 and 4 km in altitude. This corresponds to 7.5–5.5 km below the aircraft flying at an altitude of 9.5 km. Thus, the offset distance would be approximately 9.3–6.8 km, which is ≈ 0.08–0.06° of latitude. In other words, ISMAR would see 0.08–0.06° of latitude ahead of the radar (at a lower latitude) if the instruments were mounted on the same aircraft. In the simulations at 243 GHz presented here, we use the midpoint and apply a correction of 0.07°.

5.1 Results

The synthetic time series of brightness temperatures is shown in Fig. 7, along with the ISMAR data. The simulations were performed using the original aggregate models, but single oriented dendrites are also considered in Sect. 5.2. Note that the ISMAR data have been averaged over the same 200 latitude bins to get a fairer comparison (which accounts for the differences between the ISMAR measurements shown in Figs. 7c and 2d). Some of the bins contain no data, because the viewing angle is not within the relevant range of between 51 and 52°, so it was removed at the outset. This means there are regions where the lines do not join up in Fig. 7. Unfortunately, this is around the region where the largest V–H signal was measured by ISMAR.

Figure 7 Time series simulations at 243 GHz (shown in black) using a threshold value of 1700 m⁻¹ for the slope parameter λ in the exponential PSDs. Panels (a) and (b) show the H and V brightness temperatures, while panels (c) and (d) show V–H plotted as a function of latitude and V, respectively.

The simulations generally reproduce the measured brightness temperatures and V–H signal very well, in terms of both the overall magnitude and the variability. The average value of H from the measurements is 248.9 K, while the average simulated value is slightly higher at 249.3 K (Fig. 7a). The average value of measured and simulated V is the same, at 252.2 K (Fig. 7b). There is a strong correlation between the simulated and measured V–H time series in Fig. 7c, with a correlation coefficient of 0.79. Moreover, the plot of V vs. V–H in Fig. 7d shows that many of the simulations overlap the measurements. The median measured V–H has a value of 3.15 K, while the median simulated value is only slightly smaller at 2.46 K, and the interquartile range of the simulations is less than what was measured (2.82 K compared to 3.70 K). The maximum simulated V–H value of 6 K is also marginally lower than what was measured. However, considering the measurement uncertainties shown in Fig. 2d, along with our estimates of thermal noise to have an impact on V–H of about 0.5 K at 243 GHz and 3 K at 664 GHz, we do not view this as a significant error. Moreover, ISMAR and the radars are sampling a given latitude at different times. We acknowledge that sampling a given latitude at different times with the two instruments introduces error due to cloud drift. The aircraft are travelling from north to south, while the dropsonde data show that the wind direction is easterly between 1.8 and 3 km in altitude and east-south-easterly above 3 km. To calculate an estimated offset, the time difference between the ISMAR and radar measurements at a given latitude is multiplied by an estimated wind speed of 20 m s⁻¹ (Fig. 8). This shows that if one instrument samples a segment of cloud at a given latitude, then when the other instrument measures at that latitude, the feature could have moved by up to 6 km (i.e. 0.05°). In other words, after ISMAR views a certain segment of inhomogeneous cloud, the easterly winds cause the cloud to drift to the west before the radar measures at that latitude, and the cloud measured by ISMAR is not going to be seen by the radar. This explains why the results in Fig. 7 are more accurate around 58–59° latitude (where the estimated offset is within ≈ 2 km), while at latitudes where the offset exceeds 2 km, there are clear inconsistencies between the measured and simulated values.

Figure 8 The left axis shows the time difference in seconds between sampling a given latitude with the ISMAR radiometer and the 35 GHz radar. The right axis gives the estimated offset calculated by multiplying the time difference by an estimated wind speed of 20 m s⁻¹.

The sensitivity of the polarisation difference to rain rate was tested by reducing the precipitation rate from the original value of 1.12 mm h⁻¹ ($\mathit{\lambda }=\mathrm{4}\times {\mathrm{10}}^{\mathrm{3}}$ m⁻¹) to 0.5 mm h⁻¹ ($\mathit{\lambda }=\mathrm{4.7424}\times {\mathrm{10}}^{\mathrm{3}}$ m⁻¹). This was found to have very little effect on the overall results at 243 GHz increasing V–H by 0.1–0.2 K. This also implies that the melting layer is likely to be a relatively minor contributor. It is worth noting that no difference was made to the 664 GHz results by changing the rain rate. We are not getting significant emission from the lower atmosphere at 664 GHz, and results are more sensitive to mid- and high-cloud regions.

Other reasons why the largest V–H measurements are not simulated could be that we do not have the right variability in scattering as a function of size. As discussed in Sect. 4.1, an increase in simulated V–H may be possible if the dendritic monomers were oriented before aggregation or pivoted on attachment, or if the aspect ratio of the monomers was smaller. Alternatively, inclusion of single oriented dendrites may be required, which is discussed further in Sect. 5.2.

5.2 Including single oriented dendrites in the ISMAR time series simulation

The underestimation of the maximum measured V–H substantiates the theory that we need single oriented particles to simulate the strongest polarimetric signal obtained from ISMAR. In order to test that, we repeat the experiment performed in Sect. 4.2 for the time series simulations, where aggregates at the cloud base are replaced with horizontally aligned dendrites. Again we use a monodispersive distribution of dendrites with $D_{max}=$ 1 mm, with the number concentration chosen in such a way that the measured IWC in the layer is maintained. The single dendrites have a smaller aspect ratio than the aggregates, at approximately 0.1. We explore what happens when 10 %–50 % of the bottom layer of cloud is replaced. The results are shown in Fig. 9. This shows that replacing even a small proportion of the bottom layer (10 %) with horizontally aligned dendrites increases V–H, allowing the larger V–H values measured in deep cloud regions to be simulated. However, it is also clear that at some latitudes where measured V–H is lower, the simulations are more accurate without dendrites, and V–H is overestimated when they are included. Thus, large values of V–H measured from ISMAR could be an indication that there are oriented dendrites in the cloud.

Figure 9 Effect of replacing particles in the bottom layer of our model cloud with horizontally aligned dendrites. As before, red markers show ISMAR data and black markers show simulations without dendrites. The green, blue, and cyan markers show simulations when the particles in the bottom 100, 300, and 500 m of cloud are replaced with dendrites. All simulations use a threshold value of 1700 m⁻¹ for the slope parameter λ. The ISMAR data have been shifted by 0.07° latitude in panel (a) to account for the forward-viewing angle.

6 Simulating a time series of ISMAR measurements at 664 GHz

The second dual-polarisation channel on the ISMAR radiometer is at 664 GHz. Here we have performed the same time series simulations as for 243 GHz in Sect. 5 but at the higher frequency of 664 GHz. The results are shown in Fig. 10. As mentioned above, the ISMAR data were already averaged over 200 latitude bins. However, at this frequency the ISMAR data were still very noisy, so in order to reduce noise we have done a moving average over five latitude bins which corresponds to approximately 0.07° or 7.8 km. At 664 GHz, the simulations are more sensitive to the cloud layers between 5 and 8 km in altitude. Calculating the position correction in the same way as we did for 243 GHz, we apply a correction of 0.03° to the ISMAR latitude to account for the fact that the instrument is forward viewing.

Figure 10 Simulated brightness temperatures at (a) H and (b) V at a frequency of 664 GHz. Panels (c) and (d) show the brightness temperature difference V–H plotted against latitude and V, respectively. In addition to averaging the ISMAR data in 200 latitude bins, a moving average over five latitudes was applied in order to reduce noise. The same averaging was applied to the simulations for consistency. A latitude shift of −0.03° was applied to the ISMAR data to correct for the different viewing angles of the instruments.

There are clear similarities between the measurements and simulations of H and V at 664 GHz shown in Fig. 10a and b, with increasing and decreasing brightness temperatures generally predicted at the correct latitudes. The average measured and simulated value of H is almost the same (221.7 and 221.6 K). The average measured V is 222.2 K, while the average simulated V is larger at 224.7 K. There is a moderate correlation between the simulated and measured V–H time series in Fig. 10c, with a correlation coefficient of 0.67. However, the plot of V vs. V–H in Fig. 10d shows that the simulations do not overlap the measurements. There is very little polarisation signal in the measurements at 664 GHz, while the model predicts larger values up to 7.5 K. (Note that there are flights with higher 664 GHz V–H, as shown in Fig. 5 of Kaur et al., 2022.) The median measured V–H for this case has a small value of 0.35 K, while the median simulated value is bigger at 3.24 K, and the interquartile range of the simulations is greater than what was measured (4.67 K compared to 1.68 K). In terms of the separate cloud regions, the simulated brightness temperatures are quite accurate at both H and V in region a, while in region b, V is generally too large. This is particularly obvious at the start and end of this region, which is where there is a small amount of mid-level cloud visible in the 35 GHz radar (Fig. 2a). In region c, the large H depressions in regions of thicker cloud are underestimated at some latitudes and overestimated at others, while V is mainly underestimated and V–H is overestimated. The overestimation of V–H comes from a combination of errors in H and V simulations. The layer-averaged simulations in Fig. 5 show that the bottom three layers of cloud have very little effect on the H and V simulations at 664 GHz. Thus, the issues here are most likely caused by inaccuracies in the particle models in the cloud layers at 5–8 km in altitude (i.e. layers 2–4). Recall that columnar aggregate models are used in layers 2 and 3, while dendritic models are used in layer 4.

Following the discussions in Sects. 4.1 and 5.1, an underestimation in the V depression and overestimation in simulated V–H may be related to the aspect ratio of the monomers within the aggregates in these layers. The aspect ratio of the columnar monomers comprising the aggregates in layers 2 and 3 may be too big (i.e. they are too long and thin and should be more block-like), or the aspect ratio of the dendritic monomers in layer 4 may be too small (i.e. they should be thicker). It may also be that the simulations would match the measurements more closely if the monomers comprising the aggregates in these layers were rimed.

Alternatively, the small (but non-zero) measured polarisation signal suggests that there may be some element of quasi-random orientation in these layers rather than the horizontal orientation assumed in our model. This would be consistent with previous studies using CALIPSO data which have shown that oriented crystals are infrequent in cold, high clouds (e.g. Noel and Chepfer, 2010; Zhou et al., 2012). This would also agree with the findings of Fox (2020), who simulated brightness temperatures for this flight and showed that the ISMAR observations were well within the range of simulated values at 664 GHz for different ice crystal models assuming random orientation.

Another possible source of error is that we have only included particles greater than 100 µm in our model. While brightness temperature differences are not sensitive to smaller particles at 243 GHz, they are more sensitive to such particles at 664 GHz. In future work, it would be interesting to explore the sensitivity of the simulations to the degree of particle orientation and the inclusion of smaller particles in the model.

7 Conclusion

In this paper, the IMA scattering approximation was applied to a case study involving aircraft-based in situ and remote sensing observations. In situ measurements from cloud probes were used to construct a model atmosphere. Aggregates were generated, which we believe to be representative of the atmospheric conditions close to the time of the measurements. However, it is important to note that the microwave closure experiment is imperfect in design, especially for heterogeneous scenes like the one we are examining. In particular, the remote sensing and in situ measurements are not obtained at the same time, and there may be a change in microphysics between the times. There are also uncertainties arising from the different instrument viewing angles, with ISMAR viewing at approximately 50°, while the radars view at nadir. Furthermore, using in situ data has the distinct disadvantage that the cloud is only sampled in a small region, so particles may not be representative of the total cloud. Measurements are also limited to capabilities of the particular instruments, e.g. the limited sizes that can be measured by different probes. We also note that the availability of more detailed imagery would be beneficial to better constrain particle shapes.

The IMA method was used to perform scattering calculations of the generated particles. The calculations, along with the atmospheric model, were input into ARTS to perform polarised radiative transfer simulations. Comparisons of the simulated results with remote sensing measurements from the ISMAR radiometer were performed. It was found that IMA is capable of reproducing the brightness temperature depression and polarisation signature. Fox (2020) showed that simulated brightness temperatures are strongly sensitive to assumed particle shape. Here we find that choosing particle shape based on imagery generally allows for accurate simulations. The original choice of aggregates only did not fully represent the observations. It was required that some aggregates at the cloud base were changed to horizontally aligned dendrites (which were also seen in the imagery) in order to increase the V–H polarised brightness temperature differences.

In order to simulate a larger range of measurements, a synthetic time series of in situ PSDs was generated. This was done by adjusting the value of λ in the fitted exponential PSDs such that a simulation of Z matched the values measured from the 35 GHz radar. It was found that simulations of the brightness temperature depression and polarimetric V–H differences generally match the ISMAR measurements. However, the maximum simulated value of V–H is not as large as the maximum measured value. We discuss various potential reasons for this and explore some possibilities through further simulations. These include changing the intensity of rain in the model, along with adding horizontally aligned dendrites to the cloud base. We find that the polarisation difference is very sensitive to the assumed particle shape for a given ice water path, specifically the presence of single crystals mixed with aggregates. Including a small proportion (10%) of dendrites increases V–H to a more realistic value overall, while a larger proportion is required in order to simulate the largest V–H values measured. Thus, it is possible that these large polarimetric signals cannot be simulated using aggregates alone. Therefore, to obtain good retrievals from ICI, it is important to represent the split between single crystals and aggregates within the cloud as accurately as possible. Utilising the multi-frequency polarisation information available from the instrument could provide a way to constrain this, thereby reducing the need to make unrealistic assumptions.

Aside from the limitations of the experiment discussed above and the experiments we performed, other issues that are not addressed in this paper include the fact that we assume N₀ is constant in each of the seven layers for the synthetic time series simulations. In reality, it is likely that N₀ varies within each layer. Moreover, we do not represent the melting layer in our simulations. We do not include melting particles in this study since the high permittivity of liquid water means that for the scattering calculations, a large number of dipoles would be required to represent the rapid attenuation accurately. However, there was only a small change in results at 243 GHz when the rain rate was reduced, suggesting that neglecting the melting layer is not a significant source of error.

In conclusion, polarised ICI measurements are expected to provide information on oriented particles. For future missions, we recommend that dual-polarised non-nadir measurements are exploited to identify regions with oriented particles. Moreover, it would be beneficial to carry a cloud radar on the same aircraft as ISMAR.

Appendix A

The measurement uncertainties shown in Fig. 2c and d were estimated from the positive systematic error ϵ₊, negative systematic error ϵ₋, and random error ϵ_r, provided in the ISMAR data file. The uncertainties in the brightness temperatures at horizontal polarisation T_b,H were estimated as

$${\mathit{ϵ}}_{\mathrm{H}}=\sqrt{{\left({\displaystyle \frac{{\mathit{ϵ}}_{+}-{\mathit{ϵ}}_{-}}{\mathrm{2}}}\right)}^{\mathrm{2}}+{\mathit{ϵ}}_{\mathrm{r}}^{\mathrm{2}}}.$$

The region between $T_{\mathrm{b},\mathrm{H}}+{\mathit{ϵ}}_{\mathrm{H}}$ and $T_{\mathrm{b},\mathrm{H}}-{\mathit{ϵ}}_{\mathrm{H}}$ in Fig. 2c is shaded. The equivalent uncertainties were calculated for the brightness temperatures at vertical polarisation T_b,V. The uncertainties in V–H were then estimated as

$${\mathit{ϵ}}_{\mathrm{V}\text{–}\mathrm{H}}=\sqrt{{\mathit{ϵ}}_{\mathrm{H}}^{\mathrm{2}}+{\mathit{ϵ}}_{\mathrm{V}}^{\mathrm{2}}},$$

and the region between V–H +ϵ_V–H and V–H −ϵ_V–H in Fig. 2d is shaded.