Three Dimensional Radiative Effects in Passive Millimeter/Sub-Millimeter All-sky Observations

This study was conducted to quantify the errors prompted by neglecting three-dimensional (3D) effects, i.e., beam-filling and horizontal photon transport effects, at millimeter/sub-millimeter wavelengths. This paper gives an overview of the 3D effects that impact ice cloud retrievals of both current and proposed (Ice Cloud Imager) satellite instruments operating at frequencies of ≈186.3 and ≈668 GHz. The 3D synthetic scenes were generated from two-dimensional (2D) CloudSat (Cloud Satellite) observations over the tropics and mid-latitudes using a stochastic approach. By means of the Atmospheric Radiative Transfer Simulator (ARTS), three radiative transfer simulations were carried out: one 3D, one independent beam approximation (IBA), and one-dimensional (1D). The comparison between the 3D and IBA simulations revealed a small horizontal photon transport effect, with IBA simulations introducing mostly random errors and a slight overestimation (below 1 K). However, performing 1D radiative transfer simulations results in a significant beam-filling effect that increases primarily with frequency, and secondly, with footprint size. For a sensor footprint size of 15 km, the errors induced by neglecting domain heterogeneities yield root mean square errors of up to ≈4 K and ≈13 K at 186.3 GHz and 668 GHz, respectively. However, an instrument operating at the same frequencies, but with a much smaller footprint size, i.e., 6 km, is subject to smaller uncertainties, with a root mean square error of ≈2 K at 186.3 GHz and ≈7.1 K at 668 GHz. When designing future satellite instruments, this effect of footprint size on modeling uncertainties should be considered in the overall error budget. The smallest possible footprint size should be a priority for future sub-millimeter observations in light of these results.


Introduction
One of the largest sources of uncertainties in numerical weather prediction (NWP) and climate models concerns clouds containing ice [1]. The properties of such clouds are only poorly constrained in models [2] and an unmet need for trustworthy global observations of ice cloud properties exists [3].
Nowadays, satellite microwave (MW) observations are gaining weight in both weather and climate applications [4]. In fact, at the European Centre for Medium-Range Weather Forecasts (ECMWF), the assimilation of MW radiances in the so-called all-sky conditions (i.e., clear, cloudy, and precipitating) comprises 40% of the observational impact [5]. Likewise, efforts are under way to consider all-sky assimilation in high resolution regional NWP models [4]. So far, satellite instruments are supplied with MW frequencies up to ≈190 GHz. Frequencies below ≈90 GHz can penetrate nearly the entire atmosphere, and thus provide information about water vapor over ocean in clear sky conditions. In cloudy conditions, these frequencies are also sensitive to hydrometeors, with the lower frequencies (up to ≈30 GHz) giving an insight to heavy precipitation, while frequencies between 30 and 90 GHz render information about liquid water clouds. Higher frequencies between bias can be as large as 10 K. Overall, the error increases with frequency up to 37 GHz due to the increasing degree of the non-linear relationship between the brightness temperature and rainfall rate, and then, decreases at 85 GHz owing to the fact that the linear regime is reached. With the growth of vector radiative transfer models in the MW spectrum, 3D effects have also been studied in terms of polarization [35,42,43].
In all the aforementioned studies, the issue under scrutiny was the 3D radiative effects at frequencies up to 183 GHz. Only Davis et al. [44] further extended the analysis towards the sub-mm wavelengths. Based on a realistic synthetic 3D mid-latitude cirrus scene generated by a two-dimensional (2D) radar observations using a stochastic approach [45], they simulated observations from three passive instruments; namely, the Advanced Microwave Sounding Unit-B (AMSU-B) at a frequency of 190.31 GHz [46], the Earth Observing System Microwave Limb Sounder (EOS-MLS) at a frequency of 240 GHz [47], and the CIWSIR (Cloud Ice Water Submillimeter Imaging Radiometer) at frequencies 334. 65 GHz and 664 GHz [3]. They demonstrated a small horizontal photon transport effect, and a considerable beam-filling effect that affects polarization as well. However, their study is unsubstantiated, meaning that it was not supported by a statistical significance. Thus, the goal of this work is to broaden such investigations and conduce a comprehensive evaluation of the 3D effects at the millimeter (mm) and sub-mm wavelengths.
By means of the Atmospheric Radiative Transfer Simulator (ARTS) [48,49], a sensitivity study has been conducted to quantify the errors induced by neglecting the 3D structure of clouds at mm/sub-mm wavelengths. This work provides an overview of the 3D effects that impact cloud retrievals from no longer operational, current, and upcoming satellite instruments. A special emphasis is given to ICI, whose retrievals are going to provide unique information about cloud ice, with great importance for both weather and climate applications. We focus on two areas: the tropical Pacific Ocean, dominated by deep convective systems and their induced anvil cirrus clouds; and the mid-latitude North Atlantic, characterized by low pressure systems during winter. The 3D synthetic scenes have been generated from 2D CloudSat (Cloud Satellite) observations using a stochastic approach.

Instruments
This section outlines the characteristics of satellite instruments (in alphabetical order) on the basis of which we either created our synthetic scenes or we conducted forward radiative transfer simulations.

Advanced Microwave Sounding Unit-B
The Advanced Microwave Sounding Unit-B (AMSU-B) is an across-track MW radiometer comprising a channel at 89.0 GHz, a channel at 150.0 GHz, and three channels around 183.3 GHz (water vapor line) [50]. All channels are characterized by an instantaneous field of view of about 15 km at nadir. Although AMSU-B was primarily designed for water vapor sounding, the high sensitivity of its two highest channels on large hydrometeors [6,7] promoted the use of AMSU-B for precipitation [51,52] and snowfall retrievals [53]. Forward radiative transfer simulations at nadir were carried out at 186.3 GHz with a footprint size of about 15 km.

Advanced Technology Microwave Sounder
The Advanced Technology Microwave Sounder (ATMS) is a next generation across-track MW sounder that flies on board the Suomi National Polar-orbiting Partnership (Suomi NPP) in a sun-synchronous orbit at an altitude of 824 km. This instrument consists of 22 channels at frequencies between ≈23 GHz and ≈183 GHz in either vertical or horizontal quasi-polarization, with temperature, humidity, and precipitation (liquid and frozen) capabilities [54,55]. The spatial resolution at sub-satellite point is ≈75 km at 23-32 GHz frequencies, ≈32 km at 50-90 GHz frequencies, and ≈16 km at the highest frequencies. We did radiative transfer simulations at nadir at a frequency of 186.3 GHz with a 15 km footprint resembling the channel number 20 of ATMS.

Cloud Satellite
In April 2006, the sun-synchronous polar orbit satellite CloudSat (Cloud Satellite) was launched, being part of the so-called A-train constellation of five Earth observation satellites [6,56]. The payload of CloudSat is a nadir-looking cloud profiling radar (CPR) operating at 94.05 GHz that measures vertical profiles of clouds and precipitation. The vertical resolution of CPR is 480 m, but the measured backscatter signal is oversampled, leading to 125 vertical bins of about 240 m each. Radar profiles are generated every 16 ms, producing a footprint with along-and across-track resolutions of about 1.8 km and 1.4 km, respectively. In this study, the cloud geometrical profile product has been employed that includes the radar reflectivity and the cloud detection mask [6,56].

Global Precipitation Mission Microwave Imager
The GPM (Global Precipitation Mission) microwave imager (GMI) is a conical MW radiometer with 13 channels ranging in frequency from 10.65 GHz to 190.31 GHz [7]. The lowest frequencies have a footprint size of about 26 km, while the higher frequencies (above 89 GHz) have the highest-resolution footprint (6 km). All frequencies are in both vertical and horizontal polarization, except for 23.8 GHz and the channels around 183.3 GHz. GPM is a non-sun-synchronous satellite with an earth incidence angle of 52.8 • . Forward simulations are performed at 186.3 GHz with a footprint size of about 6 km.

Ice Cloud Imager
The Ice Cloud Imager (ICI) is a conical scanning radiometer that will be one of the key instruments on board Metop-SG. Metop-SG will be operated by EUMETSAT and is scheduled for launch in 2023. Following Metop, Metop-SG will have a sun-synchronous polar morning orbit with an earth incidence angle of 53 • , flying at about 832 km. ICI is going to be the first operational instrument embarking on a satellite measuring the energy emitted at mm to sub-mm wavelengths. ICI consists of 13 channels ranging in frequency from 183.31 GHz to 668 GHz, with channels ≈243.2 GHz and ≈664.4 GHz measuring both vertical and horizontal polarization, while the rest of the channels are limited to receiving vertical polarization. For all the channels, the mean footprint size, between acrossand along-track resolutions, is about 15 km. The primary objective of ICI is to monitor ice hydrometeors and humidity, filling the gap in the sub-mm wavelengths in current NWP models. Forward simulations have been conducted at frequencies of 186.3 GHz and 668 GHz at 53 • , resembling the ICI channels 183.31 ± 3.4 and 664 ± 4.2 GHz, respectively.

Microwave Humidity Sounder
The Microwave Humidity Sounder (MHS) onboard Metop-A is an across-track, five-channel MW radiometer operating at frequencies between 89 and 190 GHz, with an instantaneous field of view of 15 km at nadir. All channels are at vertical polarization, apart from the two channels at ≈183 GHz, which are at horizontal polarization. Measurements from MHS were employed originally for retrieving water vapor under most weather conditions and a few cloud products (liquid and frozen precipitation; ice water path) [52,55]. Forward simulations were conducted, as with AMSU-B and ATMS.

Special Sensor Microwave Imager/Sounder
The Special Sensor Microwave Imager/Sounder (SSMIS) is aboard the Defense Meteorological Satellite Program F-16 spacecraft, which was launched in 2003, representing the first out of a series of similar instruments that were about to be launched. SSMIS is a conically scanning passive MW radiometer with an earth incidence angle of 53.1 • which includes 24 channels at frequencies between ≈19 and ≈190 GHz. Out of the 24 channels, three channels (19.35, 37.0, and ≈91.65 GHz) include both vertical and horizontal polarization, eight channels between 57.29 and ≈63.28 GHz are right-circularly polarized, a channel at ≈22.23 GHz is only vertically polarized, and the rest of the channels are in horizontal polarization. SSMIS is a multi-purpose imager with humidity and temperature sounding channels toward an improved precipitation estimate [57]. Last but not least, the high-frequency channels enable the observation of cloud ice at a footprint size of about 15 km [58]. Hence, SSMIS has similar characteristics to the first channel of ICI, and accordingly, radiative transfer simulations are conducted at 186.3 GHz.

Atmospheric Radiative Transfer Simulator
Radiative transfer simulations have been conducted with the Atmospheric Radiative Transfer Simulator (ARTS) [48,49]. ARTS is a freely available software package, which has been designed for radiative transfer simulations for the mm and sub-mm spectral range, including polarization capabilities. ARTS has been comprehensively tested by comparison to a range of other MW radiative transfer models [59,60] and its performance has been demonstrated in various applications [11,12,14,44]. It can be operated in spherical geometry for 1D, 2D, and 3D atmospheres, and thus, it can be employed for a wide range of remote sensing applications. The key feature of ARTS concerns the availability of several scattering modules. In this study, two modules were utilized: the Discrete Ordinates Radiative Transfer Program for a Multi-Layered Plane-Parallel Medium (DISORT) [61] for 1D simulations and the Monte Carlo (ARTS-MC) [62] for 3D simulations.

Synthetic Scenes
The synthetic 3D cloud fields were generated by an Iterative Amplitude Adjusted Fourier Transform (IAAFT) algorithm [63]. This algorithm is applied on 2D CloudSat measurements (pressure, latitude) of radar reflectivity (in dBZ) and yields so-called surrogate 3D fields (pressure, latitude, and longitude) that share the amplitude distributions and the power spectra of the measured fields. This algorithm works by iteratively adapting the amplitude distribution and the Fourier coefficients to obtain the desired power spectra. For details, the reader is referred to [63,64]. From the resulting 3D fields of radar reflectivity, the ice water content (IWC) and rain water content (RWC) are retrieved (see Section 2.5). This algorithm was previously utilized and its performance was successfully tested [65]. In brief, Rydberg et al. [65] developed a Bayesian retrieval algorithm for obtaining cloud ice and water vapor from Odin-SMR limb measurements, with CloudSat measurements being the basis for deriving 3D cloud structures.
The remaining required atmospheric and surface components, i.e., water vapor, temperature, liquid water content (LWC), and surface temperature, were extracted from collocated global atmospheric reanalysis data produced by ECMWF, the ECMWF reanalysis data (ERA-Interim). Note here that water vapor and liquid water content profiles were scrutinized. In order to avoid unrealistic values of water vapor in areas with clouds or precipitation, a value of 0.9 (in terms of relative humidity) was employed in areas with values of hydrometeor content (i.e., RWC or IWC) above 50 µg m −3 . Moreover, LWC was set to zero at temperatures lower than 243 K and at heights above 10 km. Last but not least, in bins with no retrieved hydrometeors, liquid water content was set to zero, ensuring there was no artificial placement of liquid clouds from ERA-Interim that were not detected by CloudSat.
In order to conduct a comprehensive statistical analysis and quantify the magnitudes of both the beam-filling and horizontal photon transport effects, we randomly sampled CloudSat scenes. In the tropics, CloudSat overpasses within N5 • to N9.4 • for July 2015 were selected. Similarly, in the mid-latitudes, CloudSat overpasses within N48.8 • to N53.2 • for January 2015 were chosen. This selection led to a total number of 30 CloudSat overpasses in the tropics and 29 overpasses in the mid-latitudes. Then, each orbit was further divided by 2.2 • leading to 55 (58) scenes over the tropics (mid-latitudes). By performing IAAFT, the resulting 3D scenes were characterized by a grid size of about 160 km by 200 km. Figure 1a illustrates an example scene used as a basis for the generation of the synthetic 3D field. It corresponds to the 2D radar reflectivity field measured by CloudSat for an overpass within N5 • to N7.2 • , where a deep convective cloud system with an overlying thin cirrus layer and a small convective core were identified. An example latitudinal slice of the synthetic radar reflectivity field and the maximum of the surrogate radar reflectivity for each latitude and longitude generated by IAAFT were found in Figure 1b and 1c, respectively. The highlighted black box denotes the area where the actual forward radiative transfer simulations were conducted (see Section 3).

Hydrometeor Microphysical Properties
Accurate radiative transfer simulations of brightness temperatures necessitate realistic representations of the hydrometeor's size, density, shape, and size distribution. Although several studies focused on the impacts of the microphysical properties on brightness temperatures simulated at low frequencies (up to ≈200 GHz) [66][67][68][69][70], limited research has been conducted in sub-mm wavelengths [71]. Even though such studies are of great importance on their own, they are beyond the scope of this paper.

Particle Size Distributions
Two particle size distributions (PSD) have been considered in this study, describing the distributions of the different particle sizes per volume element. In the tropical region, the PSD by McFarquhar and Heymsfield [72], hereafter MH97, was utilized, which is determined by IWC and temperature. MH97 induces higher loads on small particles and has been proven a realistic PSD for tropical conditions [47,73,74]. In the mid-latitude region, the PSD developed by Field et al. [75], denoted as F07, was employed. F07 is specified by temperature, IWC, and the mass-size relationship and has been demonstrated for use in mid-latitude conditions by in situ observations [12,68,76]. In comparison to MH97, F07 exhibits a stronger preference towards larger hydrometeors with increasing ice water content. Last but not least, the PSD by Wang et al. [77] was used to describe the rain water content.

Single Scattering Properties
Selecting the optical properties of ice hydrometeors for radiative transfer simulations is quite troublesome due to their high heterogeneity in size, shape, and orientation within a cloud system. In this work, the ARTS scattering database has been considered [78]; it offers the single scattering properties of more than 30 hydrometeor shapes, frequencies (1 to 886.4 GHz), and sizes, and thus, meets our requirements.
In the framework of this study, different ice hydrometeors' shapes have been tested, but their sensitivity to the 3D radiative effects was found to be of low importance. Consequently, a comprehensive assessment of the latter effects have been conducted only for large plate aggregates that are found to be one of the realistic scatterers at sub-mm wavelengths [71]. For liquid particles (cloud water and rain), the optical properties have been derived from Mie theory in the case of spherical particles.
Note here that the orientation of hydrometeors is also of great importance, but due to the lack of scattering databases that cover orientation, only randomly oriented hydrometeors have been considered.

Bulk Scattering Properties
The bulk scattering properties are derived from the single scattering properties by averaging over the hydrometeor number density (HND). HND is calculated by the hydrometeor content (HC) and temperature information, i.e., IWC (below 273.15 K) and RWC (above 273.15 K), by integrating the assumed PSD over all sizes. In other words, it describes the number of hydrometeors per size bin.

Radar Inversions
For a given particle size distribution N(D eq ), with D eq describing the volume equivalent diameter, and a hydrometeor shape, the radar reflectivity factor Z e (in units of mm 6 m −3 ) is given by [79]: where c is the speed of light, ν is the frequency, K denotes the dielectric factor [56], and σ is the back-scattering cross section. Z e is usually expressed in logarithmic scale: Z e (dBZ) = 10 · log 10 Z e (mm 6 m −3 ).
Similarly, the hydrometeor content is defined in terms of the particle size distribution and the mass-size relation m [80]: From Equations (1) through (3), the relationship between the hydrometeor content and radar reflectivity is inferred; i.e., HC-Z e , and the HC are retrieved for a given microphysical combination (PSD and hydrometeor shape). Following an iterative approach (onion-peeling) [81], the HC is estimated from the equivalent Z e from each atmospheric layer, starting from the uppermost to the lowermost one. At each iterative step, the retrieved value of HC is corrected to account for attenuation by hydrometeors and atmospheric gases. Note here that IWC and RWC are estimated separately, with a temperature of 273.15 K being the decisive indicator. For details, the reader is referred to [71]. Figure 2 illustrates the IWC field corresponding to the example latitudinal slice of the synthetic surrogate radar reflectivity field depicted in Figure

Radiative Transfer Simulations
By means of ARTS, forward pencil beam (infinite width) radiative transfer simulations have been carried out at two frequencies, i.e., 186.3 GHz and 668 GHz; two incident angles at the ground, i.e., 53 • (slant view) and 0 • (nadir view); and two footprint sizes, i.e., 6 km and 15 km, representing the characteristics of no longer operational (AMSU-B and SSMIS), current (AMSU-B, ATMS, GMI, and MHS), and proposed (ICI) satellite instruments. Note here that for all the across-track viewing instruments (AMSU-B, ATMS, and MHS), the nominal footprint size at nadir is considered and any mention to nadir view should not cause any confusion. Although a microwave radiometer operating at a frequency of 668 GHz and characterized by a footprint size of 6 km does not exist yet, we decided to include it into the analysis as a basis for comparison to the aforementioned instruments and as an indication for future satellite missions. In this study, simulations were conducted at only two frequencies covering the margins of the frequency range at which ice cloud retrievals are or are going to be based on. That way, we could draw conclusions for the magnitude range of the 3D effects. Covering all the frequency channels of ICI was beyond the scope of this paper and requires a dedicated study.
The absorption coefficients and vertical profiles of the major atmospheric gases have been considered as follows: the gas absorption model PWR-98 [82] was used for water vapor and oxygen, the model by Liebe et al. [83] was selected for nitrogen, and the liquid absorption model of Ellison [84] was utilized for liquid water clouds. All radiative transfer simulations were conducted over ocean. In that direction, the recently developed parameterization of sea surface emissivity generated by the Tool to Estimate Sea Surface Emissivity from Microwave to sub-Millimeter waves (TESSEM2) [85] was selected. Moreover, lambertian reflection was assumed. In fact, the method of handling the surface contribution to the simulated brightness temperature is of little concern since, for the frequencies under investigation, it is negligible [10].
For each scene, three radiative transfer calculation modes were used: 1. A realistic full 3D pencil beam simulation (3D mode) using ARTS-MC at a 2 km spatial grid, whereby a target precision of 1 K was selected (MC noise). Note here that in order to avoid excessively long simulations linked to domain areas with enhanced multiple-scattering events, the maximum number of allowed photons was set to 10 5 . 2. An independent beam approximation (IBA mode) denoting the independent column approximation (ICA) employing the DISORT solver. Accordingly, 1D radiative transfer simulations were conducted along the line of sight at a 2 km spatial grid, accounting for horizontal heterogeneities, but omitting any horizontal photon transport (HPT). 3. IBA simulations following the plane-parallel assumption; i.e., utilizing domain-averaged bulk scattering properties (1D mode). This averaging was conducted by applying a 2D Gaussian function (antenna pattern) characterized by full widths at half maxima of 6 km and 15 km, which represent the two footprint sizes of the aforementioned satellite instruments. For each 1D mode, two individual simulations were performed according to which step we applied the antenna pattern of the sensor: (a) HND-avg, where we directly averaged the 3D HND field and (b) HC-avg, where we first averaged the HC, and then derived HND as described in Section 2.4.3.
For all the calculation modes, the simulated brightness temperature (T B ) is expressed in terms of the hydrometeor impact: where T B,clear denotes the brightness temperature for the same atmosphere, but in the absence of clouds. Figure 3 illustrates the simulated ∆T B for pencil beam calculations at 668 GHz for the full 3D mode Figure 3a, IBA mode Figure 3b, and the differences between the latter calculation modes Figure  3c. This is an overview of the simulated brightness temperatures at the original 2 km grid. A rather good agreement was found between the two calculation modes, with the IBA mode leading to slightly smoother results compared to the 3D mode. However, in areas with large gradient in hydrometeor content (see Figure 3c), the simulated differences between the two modes were up to ± 15 K. For the 3D and IBA pencil beam calculations at a 2 km spatial grid, the antenna pattern of the different sensors was applied post-simulation. In other words, we integrated the simulated brightness temperature over the field of view of the sensor considering 3D heterogeneous atmosphere, but neglecting HPT in the case of IBA mode. Henceforth, any mention to 3D or IBA mode corresponds to area averaged results (over the antenna pattern). Following Marshak et al. [86], we considered the following effects: By comparing 3D and IBA mode, the significance of the horizontal photon transport effect was assessed: Confronting IBA with 1D, we rendered information about the importance of the horizontal heterogeneity effect; i.e., beam-filling effect (BF): Accordingly, the 3D simulated brightness temperature minus the 1D counter part results leads to the total 3D effect: A first indication of the HPT effect is shown in Figure 3c, but by means of pencil beam simulations (without considering the antenna pattern). Note here that by applying the 2D Gaussian function to the pencil beam calculations, the differences between the two modes will be considerably smoothed out. In addition, the MC noise (1 K) is significantly decreased; i.e., at least by an order of magnitude.   Figure 1c: (a) generated by the 3D mode, (b) generated by the independent beam approximation (IBA mode), and (c) differences in the simulated brightness temperature between 3D and IBA.
As it was mentioned is Section 2.4.2, only randomly oriented hydrometeors are considered, and thus, the issue under scrutiny is the 3D radiative effects in the scalar radiative transfer theory, i.e., only the first Stokes component is simulated and polarization is neglected.

Results
Results are presented for all synthetic scenes, which were based on CloudSat overpasses over the tropics for July 2015 and over mid-latitudes for January 2015. Footprints with a hydrometeor impact less than −3 K were excluded from figures and statistics; varying this limit by ≈1 K leads to differences below 0.5 K in terms of the root mean square error (RMSE). As a result, the size of the sample for each frequency and footprint size varies, with the smaller sample being related, primarily, to the lower frequency. Tables 1-3 summarize the statistics for the horizontal photon, the beam-filling, and the total 3D effect, respectively, while Figures 4-7 highlight the findings in the case of the intertropical scenes, wherein the largest 3D effects are found. To highlight the differences in the magnitude of the BF effect between tropics and mid-latitudes, Figure 8 illustrates an example.

Tropics
An illustration of the direct comparison between 3D and IBA simulations for a slant view of the different sensors is shown in Figure 4. This comparison reveals only small scatter, with a bulk of data points aligned along the one-to-one; the one-to-one line (black line) and the linear fit (dark green line) lie almost on top of each other, with the slope of the linear fit (α) being 0.988 < α < 0.991. Only at large values of hydrometeor impact (≈ ∆T B < −75 K) can one slightly tell apart the two lines, especially at the frequency of 668 GHz for both footprints. Overall, an excellent linear correlation was found between the two modes, with the larger scatter being evident, for a given sensor footprint size, at a frequency of 668 GHz compared to 186.3 GHz, and for a given frequency, at the smallest footprint size. The RMSE was below 1.964 K in all cases. Note here that the number of photons employed for the 3D mode was sufficient to yield accurate simulated T B , leading to very small MC noise. Solely for a sensor footprint size of 6 km, the two-standard deviation (2σ) of the MC is somewhat visible.  A more comprehensive analysis is depicted in Figure 5, which corresponds to the histogram of the differences in brightness temperature between the two calculation modes for slant view. Hence, it provides an overview of the magnitude of the horizontal photon transport effect. For all the frequencies and footprints, the differences between 3D and IBA simulations are roughly symmetric, are clustered close to zero (marked by the black line), and can be described quite well by normal distributions. Only small discrepancies were found, with IBA simulations slightly overestimating the simulated T B and leading to a mean HPT error below 0.433 K. However, apparent outliers exist (see the crosses in Figure 5), which are mostly evident at the frequency of 668 GHz and can be as large as ≈9 K (99 th percentile; not shown here). From nadir, the magnitude of the HPT effect is less pronounced compared to that from the slant view. By conducting simulations at nadir, the photon path length decreases significantly, and accordingly, the number of scattering events decreases. In brief, the mean error is up to 0.226 K with a RMSE up to 0.676 K. The linear correlation between the two calculation modes is excellent (R 2 ≈ 1), indicating that almost 100% of the variance in the simulated brightness temperature derived by the full 3D simulations can be described by IBA counter part simulations. For details, the reader is referred to Table 1.

Mid-Latitudes
In mid-latitudes, the HPT effect is smaller than in the tropics, with a mean error below 0.148 K and a RMSE of about 0.330 K at both frequencies, footprints, and instrumental views.
For radiative transfer simulations from both slant and nadir views, the histograms of the HPT effect (not shown here) at 668 GHz are centered closer to the zero line compared to 186.3 GHz, leading to a smaller mean error (see Table 1). However, in the case of slant view, the corresponding RMSE was found to be larger at the highest frequency and the smallest footprint size, in agreement with the tropics. From nadir view, although the mean differences between 3D and IBA simulations were found to be larger with the smaller footprint size, the magnitude of the HPT errors are comparable among the frequencies under investigation (differences are in the second decimal point). Only for 186.3 GHz, was the mean HPT slightly larger at 15 km (but not in terms of the RMSE). Moreover, comparing the HPT effect from the slant and nadir views, it follows that for the former case, the histograms are centered closer to the zero line leading to a slightly smaller HPT error. Table 1. Statistics of the horizontal photon transport effect over both tropics and mid-latitude scenes for both slant and nadir radiative transfer simulations. Results are presented in terms of the sample size (n), mean, root mean square error (RMSE), slope of linear regression line (a), two-standard deviations of the regression slope (2σ), the coefficient of determination (R 2 ), and the RMSE corrected by the slope α of the regression line (RMSE c ). The correlations of the simulated brightness temperature between the IBA and the two different 1D calculation modes (see Section 3) at slant view, i.e., HND-avg (in dark blue) and HC-avg (dark red), are found in Figure 6. In contrast to the comparison between the 3D and IBA simulations (see Figure  4), this comparison introduces a significant scatter, especially at a frequency of 668 GHz. Still, most of the points are bundled around the one-to-one line, mainly in the case of HC-avg (and for ∆T B below ≈ −80 K) and rather good linear relations were found between the IBA and 1D, with R 2 values being greater than 0.957. Altogether, performing 1D simulations led to an overestimation of the hydrometeor impact, with the overestimation found to be larger at the highest frequency (668 GHz) and at the largest footprint size (15 km). Only for a footprint size of 6 km and a frequency of 186.3 GHz, was a rather good agreement found between the two calculation modes, where the smallest scatters were found with R 2 values of about 0.997 and 0.999 for HND-avg and HC-avg, respectively. Furthermore, between the two 1D calculations modes, discrepancies are larger for HND-avg compared to HC-avg. Note here that in the case of HC-avg, for rather low values of the hydrometeor impact (above −20 K), the data points are above the one-to-one line.  . The same as in Figure 4, but between IBA and both 1D calculation modes; i.e., hydrometeor number density average (HND-avg in dark blue) and hydrometeor content average (HC-avg in dark red) modes. Figure 7 displays the corresponding histograms of the differences in brightness temperature between IBA and 1D calculations modes (at slant view) and gives an insight into the beam-filling effect. The discrepancies associated with the HND-avg (in dark blue) are concentrated on the left, and hence, are characterized by a right-skewed distribution, with a rather long tail, for both frequencies and footprints. Additionally, the peaks of the distributions lay always to the right from the zero difference line (marked by the black line), signifying a consistently positive beam-filling effect, leading to a mean error up to 1.549 K at 186.3 K and one up to 5.966 K at 668 K; the corresponding RMSEs are 3.973 K and 13.13 K, respectively. To summarize, the beam-filling error induced by neglecting domain heterogeneities from the slant view can be as large as ≈14 K at the highest frequency, defined by the 1.5 times the interquartile range (see the right arm of the whisker box in Figure 7d). Note here that the outliers can be as large as ≈40 K at 168.3 GHz and even more than 50 K at 668 GHz, but in the graphs results are presented up to a BF effect of about 15 K for better illustrative reasons.

Tropics
In the same direction are the findings for the HC-avg at 668 GHz; i.e., a right-skewed distribution was identified-characterized by a long tail towards the right. However, the magnitude of the BF effect is less pronounced, with a mean error (RMSE) of 1.550 (5.737) K at a footprint size of 6 km and an RMSE of 4.070 (10.76) K at a footprint size of 15 km. On the contrary, at a frequency of 186.3 K the histograms of differences in T B approximate a normal distribution quite well, with a very short tail (especially at a 6 km footprint size). Accordingly, the mean (RMSE) value of the BF is −0.065 (1.047) K at 6 km footprint size and is 0.160 (2.123) K at 15 km. For both frequencies and footprint sizes in Figure 7, one can clearly see that, in this case (HC-avg), the peaks of the histograms are always located to the left of the zero difference line. In other words, at low values of the hydrometeor impact (above −20 K), conducting 1D radiative transfer simulations by means of a HC-avg approach can lead to a negative BF effect (an underestimation of the simulated brightness temperature).
For radiative transfer simulations at nadir, the same behavior was found. However, conducting 1D simulations (both HND-avg and HC-avg) led to slightly smaller BF errors compared to the simulations at slant view. From nadir, the field of view of the sensor is smaller (on average), and therefore, so is the heterogeneity along this path. Consequently, slightly better linear correlations were found between IBA and both 1D simulation modes, with 1D simulations explaining above ≈96% and 97% of the variance in the simulated ∆T B for HND-avg and HC-avg, respectively.

Mid-Latitudes
In the mid-latitudes, the errors induced by neglecting domain heterogeneities followed the same trends as in the tropics; i.e., the errors were found to be larger at the largest footprint size (15 km), at the highest frequency (668 GHz), and from slant view rather than nadir. Thus, we report only the slight differences, and for details the reader is referred to Table 2. To begin with, compared to the tropics, the simulated hydrometeor's impact is considerably lower, with its values being above −40 K and −60 K at 186.3 GHz and at 668 GHz, respectively. Clouds over mid-latitudes were found at lower altitudes, and thus, were characterized by lower ice hydrometeor content. Consequently, the beam-filling effect was found to be significantly smaller with a mean error below ≈0.9 K and a RMSE below 2.437 K. Regarding the two 1D calculation modes, both HND-avg and HC-avg produce comparable BF effects, with HND-avg leading to the largest error, only at 186.3 GHz. At 668 GHz, the opposite holds true.
The histograms of the BF effect are described by narrower distributions compared to the case of tropics. As an example, Figure 8 depicts the corresponding histograms in the case of simulations from slant view. For both footprint sizes and both frequencies, only right-skewed distributions were found for both HND-avg and HC-avg. Recall here that in the case of the tropics and a frequency of 186.3 GHz, the corresponding histograms for HC-avg were better represented by normal distributions. Considering the low hydrometeor impact and the small magnitude of the BF error, the probability density at smaller values is high compared to that in the tropics (see the peaks for each mode in Figure 7). An exception is posed by the distribution of the BF effect at a frequency of 186.3 GHz with a footprint size of 15 km in the case of HND-avg only, where, apart from the high peak that is located left from the zero difference line, a much lower peak is clearly visible at values of the BF error close to 1.6 K. The same behavior is visible in the histograms of the BF error from nadir view, but the two peaks are slightly closer to each other (not shown here). For the two 1D calculation modes, both HND-avg and HC-avg produce comparable beam-filling effects, with HND-avg leading to the largest error only at 186.3 GHz. At 668 GHz, the opposite holds true.

Total Effect
As we already mentioned, the comparison between the full 3D and the 1D simulations introduces the total effect; i.e., the combined impact of the horizontal photon transport and beam-filling effect. Table 3 compiles the corresponding statistics. One can clearly see that on average, the total effect can be seen as the sum of the two effects (comparing Table 3 with Tables 1 and 2), but this is no surprise considering their linear relation (see Equation 7). For very small values of the hydrometeor impact (below ≈ −10 K), there is a rather good agreement between the calculation modes. However, the differences increase at higher values of the ∆T B . Overall, the total effect follows the same behavior as the BF effect, and consequently, any repetition is redundant. For details, the reader is referred to Table 3 and Section 4.2. From this, in addition to the small horizontal photon effect, one can conclude that the total radiative effect is always driven by the BF one.
On the whole, conducting 1D simulations and neglecting the 3D heterogeneous atmosphere over the field of view of the sensor can lead to substantial 3D effects, with a mean error (RMSE) up to 1.738 (4.425) K at 186.3 GHz and one up to 6.337 (13.88) K at 668 GHz.  Table 3. The same as in Table 2, but in terms of the combined beam-filling and horizontal photon transport effect; i.e., the total 3D effects (Total = ∆T B (3D) − ∆T B (1D) = HPT + BF).

Empirical Correction Scheme
Each of the aforementioned comparisons among the different calculation modes, i.e., 3D against IBA (HPT effect), IBA against 1D (BF effect), and 3D in contrast to 1D (total effect), can be described by a linear regression model. To that end, the ordinary least squares regression technique was utilized, whereby the following simple model function was applied: This technique yields the projection line, whereby each of the regressors x (IBA or 1D) are correlated to a maximal degree with y (3D or IBA) and supply the slopes, which reduces the error in the prediction of y. Results are compiled in Tables 1-3. As described above, an excellent linear correlation was obtained from all the comparisons, but with a significant scatter, mainly in the case of the total and the beam-filling effects, which strongly depend on the frequency, footprint size, and slant path. To that end, the slope α of the latter model can be seen as a multiplication correction factor (hereafter slope correction) that forces the errors in the prediction of the simulated brightness temperature by each approximation (IBA or 1D) to be more symmetric. Note here that the impact of such a slope correction at mm/sub-mm wavelengths is more pronounced when the 3D effects are larger. Accordingly, at the lowest frequency, smallest footprint size, and from nadir view, it induces only a minor improvement of about 0.1 K. Starting with the total effect, both 1D simulation modes can explain at least 95.2% (see Table 3) of the variance in the brightness temperature field resulting from the full 3D mode, with the HC-avg method leading to slightly better results. Consequently, the slope correction diminishes the RMSE by up to 3.2 K, chiefly at a footprint size of 15 km and a frequency of 668 GHz, where the largest effect was found. Overall, the RMSE of the combined effect introduced by conducting 1D simulations was improved from 0.434 < RMSE < 13.88 K to 0.351 < RMSE c < 10.65 K in the case of HND-avg and from 0.386 < RMSE < 11.53 K to 0.332 < RMSE c < 1.854 K in the case of HC-avg, where RMSE c denotes the corrected RMSE. In the same direction are the findings for the beam-filling effect, in which the slope correction can lead to an improvement up to ≈3 K. On the contrary, the slope correction is negligible (below 0.1 K) when it comes from the horizontal photon transport effect. In fact, IBA simulations are able to explain nearly 100% of the variance in the 3D simulations, and only minor discrepancies were recorded (see Section 4.1). For details, the reader is referred to Tables 1-3.

Discussion
By confronting the full 3D with 1D simulations, one can explore the errors induced by neglecting the 3D aspect of the atmosphere and give rise to the so-called 3D effects-the beam-filling and the horizontal photon transport effect. In order to provide estimates of the magnitudes of these two effects individually and point out which has the dominant effect at mm/sub-mm wavelengths, we additionally conducted simulations according to the independent beam approximation.

The Horizontal Photon Transport Effect
At mm/sub-mm wavelengths, the single scattering albedo of ice hydrometeors is very high, meaning that their interaction with radiation is mostly governed by scattering. In other words, the multiple scattering process is enhanced, while absorption, and hence, emission, is of less importance. Accordingly, one could assume that the horizontal photon transport effect, which is driven by a multiple-scattering process, would be enhanced compared to other studies at lower frequencies, where the issue under investigation was liquid water clouds or precipitation that are characterized by smaller single scattering albedo [39,87]. On the contrary, the comparison between the full 3D radiative transfer simulations (3D mode) and the independent beam approximation (IBA mode) reveals no considerable discrepancies, with IBA simulations introducing a small overestimation and mostly random errors. In short, the mean HPT error is below 0.433 K, with a RMSE below 1.964 K, confirming the findings of Davis et al. [44] over the same frequencies. Particularly, IBA simulations (without considering the antenna pattern) produce somewhat smoother results compared to 3D simulations that could lead to both positive and negative horizontal photon transport effects. These effects, which can be as large as ±15 K, e.g., see Figure 3c, are subject to the 3D structures of clouds, especially the areas with a large gradient in hydrometeor content. However, the integration over the antenna pattern of the satellite, which is characterized by a large full width at half maximum (6 and 15 km), neutralizes the effect and significantly decreases its size.
The magnitude of this effect strongly depends on the frequency, the footprint size, and the observation angle of the satellite. In fact, the larger the frequency, the larger the effect. Larger frequency means that the hydrometeors are characterized by a larger scattering cross section, which potentially leads to an increasing number of multiple-scattering events, and thus, an increasing probability for a photon traveling between neighboring columns with different optical properties. In addition, decreasing the footprint size results in an increase of the HPT effect because the variability in optical properties is resolved by smaller scales. In principle, the magnitude of the HPT effect is less pronounced from nadir view compared to slant view. This can be explained by the fact that the effective photon path length is (on average) decreased, and consequently, so is the number of scattering events. However, this is true only in the tropics. In the mid-latitudes, the mean errors owing to the neglection of the HPT were found to be slightly larger for nadir viewing instruments. This can be explained by the fact that the differences between the 3D and IBA modes are centered closer to the zero difference line (see Figure 5). In terms of the RMSE, the HPT effect is still larger for the slant path-with those instruments operating at the lowest frequency and being characterized by the largest footprint size (ICI-and SSMIS-wise) being the exception; recall here that for all across-track instruments, we only refer to their individual nominal resolutions at nadir. For the solar spectrum, similar results have been reported by Barlakas [22] and OHirok and Gautier [88]. They found that such behavior is subject to the gradient in hydrometeor content and the observational geometry and could be attributed to photon trapping and shadowing effects. Last but not least, the HPT in mid-altitudes is less pronounced compared to the tropics. Clouds containing ice in mid-latitudes are located at lower altitudes and are characterized by lower ice water content. In addition, at lower altitudes, the absorption by water vapor is of great importance for the frequencies under investigation, and hence, the number of scattering events is rather low, which further diminishes the magnitude of the HPT effect.

The Beam-Filling Effect
For both tropics and mid-latitudes, the beam-filling errors induced by conducting 1D radiative transfer simulations from both slant and nadir views increase with increasing frequency. However, increasing the footprint size of the sensor results in an increase of the magnitude of the BF effect. This can be explained by the fact that we resolve the variability in optical properties at larger scales. Overall, performing 1D simulations over all-sky conditions leads to a systematic overestimation of the hydrometeor impact, with mean error (RMSE) up to 1.549 (3.973) K at 186.3 GHz and one up to 5.966 (13.13) K at 668 GHz. This overestimation, which is driven by the non-linear relationship between the brightness temperature and the hydrometeor content in the field of view of the sensor [3], further supports the findings of Davis et al. [44].
At rather low values of the hydrometeor impact, i.e., below ≈ −20 K, the errors prompted by neglecting domain heterogeneities can produce both positive and negative beam-filling effects (see Figure 6). At the frequencies under investigation, the effect's sign depends on the vertical and horizontal variability (gradient) of the hydrometeor content, and consequently, the vertical and horizontal placement of clouds along the field of view of a sensor. In brief, by averaging inhomogeneities over the large footprint sizes, especially at slant view instruments (AMSU-B-, ATMS-, and MHS-wise), could potentially lead to values of ice water path that can be high enough to be at the non-linear regime, and 1D simulations could yield higher hydrometeor impact compared to IBA simulations for a given ice water path [3]. Similar findings have been found by Davis et al. [44] and Lafont and Guillemet [36]. Davis et al. [44] reported a negative BF effect in the case of mm wavelengths (190 GHz) at low values of the hydrometeor impact (below ≈7 K). In addition, Lafont and Guillemet [36] demonstrated that for convective precipitating clouds, the BF effect can be between −25 to 10 K strongly depending on cloud cover. Between HND-and HC-avg calculation modes, the latter one mostly leads to negative values of the BF effect and to an overall smaller bias. One can clearly see in Figures 7 and 8 that the peak of the histograms for the latter case are always located to the left of the zero difference line. Averaging the hydrometeor content over the field of view of the sensor, and then determining the hydrometeor number density, could potentially induce particles that are smaller in size. Consequently, for a given frequency, this implies particles characterized by smaller scattering cross section. In other words, HC-avg results in hydrometeors that are more absorbing (enhanced emission) and by confronting them to the IBA counter part brightness temperature could lead to the negative sign of the BF effect. For instance, by directly applying the antenna pattern to the hydrometeor number density (HND-avg), we do not alter the size of hydrometeors, and therefore, the impact of the aforementioned non-linearity is more pronounced.

Total Effect and Implications for Retrievals and Data Assimilation
The total effect prompted by ignoring 3D heterogeneities was found to be consistent with the beam-filling effect. Moreover, with the paucity of the horizontal photon transport effect, we could conclude that, for mm/sub-mm wavelengths, the beam-filling one drives the 3D radiative effects.
Our analysis implies that observations of cloud ice, and as a result, ice retrievals based on all the aforementioned instruments, are and will be unfavorably contaminated by the beam-filling effect. As we mentioned above, the BF effect, and thus, the 3D effects overall depend on the view (slant or nadir), the size of the footprint, and the frequency, with the latter one causing the largest impact. Thus, ICI-based retrievals are the ones mostly affected by 3D effects characterized by RMSEs of up to ≈4.4 K at 186.3 GHz and ≈14 K at 668 GHz. The same stands for any future satellite mission with instruments operating at sub-mm wavelengths, which are described by even smaller footprint size (6 km), whereby the corresponding RMSE could be as large as ≈8.2 K. Among the current instruments that operate at a frequency of 186.3 GHz and are characterized by a footprint of about 15 km (AMSU-B, ATMS, and MHS at nadir only and SSMIS), SSMIS-based retrievals are marginally subject to the largest error (up to ≈1.2 K more). GMI-wise retrievals, which are supplied by a channel at ≈186.3 GHz and characterized by the smallest footprint size, are the ones affected least by 3D effects; the RMSE is below ≈2.6 K. These results allow us to draw similar conclusions for the upcoming Microwave Imager (MWI) mission that will fly together with ICI on Metop-SG. This instrument is a conically scanning radiometer with five channels at ≈183 GHz (among others). Even though MWI is characterized by a footprint of about 10 km, we assume that the beam-filling RMSE will be around ≈3.5 K, especially in scenes over intertropical latitudes.
In passive MW retrievals of precipitation, the errors induced owing to the use of 1D plane-parallel simulations have been investigated by plenty of former studies. Consequently, correction schemes have been developed over the last years consisting of a static factor [30,36], a dynamic adjustment [89], and diminishing the error by employing the channel with the highest resolution [34,38]. In this study, a first attempt was made to mitigate 3D effects by applying a statistical correction scheme by means of a multiplication factor. On the basis of the excellent linear correlation among the different calculation modes, we wield the ordinary least squares regression method to derive the slope of the simple linear model as the multiplication factor to induce a more symmetric 3D error. Accordingly, the slope correction diminished the RMSE by up to 3.2 K, particularly at the largest frequency and footprint size.
Finally, our findings are also relevant to all-sky data assimilation. In data assimilation, sub-grid variability (domain heterogeneities) is only crudely taken into account. Accordingly, radiative transfer simulations are conducted by means of a 1D plane-parallel approximation. In addition, considering the much larger scales of the model domain, a considerable 3D effect is expected. For instance, for the frequencies considered so far (below ≈190 GHz), the area-averaged precipitation is systematically underestimated [90]. ICI holds a strong potential to extend the scope of data assimilation towards the sub-mm region. Accordingly, our results provide a first insight of the size of the 3D effect and further explore the use of an empirical correction scheme, which could potentially yield a more symmetric error.

Implications to Future Satellite Instruments
The significant 3D effects discussed above should be considered in the design of future satellite missions. Let us consider, for example, the highest frequency channel of ICI; i.e., 664 GHz. For this channel, the noise requirement is 1.6 K [10]. For observations over tropical latitudes that are substantially affected by clouds, there is also a modeling uncertainty of ≈9.5 K, even if a correction is applied (see Table 3). The two error sources discussed are independent and the total uncertainty is ≈9.6 K.
A MW radiometer operating at the same frequency, yet with a much smaller footprint size, would be subject to smaller modeling uncertainties. For the case of a 6 km footprint, the matching estimate in Table 3 is 6.0 K. An instrument with such a footprint size can be achieved. In fact, the 664 GHz channel of ICI is designed to have the same footprint size as the 183.31 GHz ones, while the antenna size, in principle, allows for a footprint of 6 km. For an across-track scanning instrument, an effective antenna size of 1 dm would suffice to obtain 6 km resolution over a considerable swath width. Anyhow, assuming that other technical characteristics are unchanged, observations for a 6 km footprint would have a noise level of 4 K (a factor of 2.5 higher), leading to total uncertainty of 7.2 K.
Overall, a more accurate cloud retrieval can be performed at 6 km resolution than at 15 km, despite that the instrument noise is higher. One could argue that a noise of 4 K is too high for a meaningful retrieval of humidity. However, for clear-sky conditions the 6 km observations can be averaged to 15 km resolution and the noise is brought back to 1.6 K. An instrument having a 6 km resolution will require a faster rotating antenna, which will cause technical challenges, and will generate larger data volumes. Nevertheless, aiming for smaller footprint sizes should be considered, since it will provide cloud information with both higher spatial resolution and better precision.
Future satellite missions could also explore the combined use of an active instrument (e.g., similar to CloudSat) in a constellation with a passive radiometer operating at sub-mm wavelengths. The high vertical and horizontal resolution of the active sensor could supply the passive retrievals with a more detailed view of sub-grid variability of ice hydrometeors. Such a combination gives, in addition, increased information on microphysical properties of ice hydrometeors [91].

Conclusions
By means of the Atmospheric Radiative Transfer Simulator (ARTS), a study has been carried out in order to explore the errors prompted by omitting three-dimensional (3D) radiative transfer at millimeter/sub-millimeter (mm/sub-mm) wavelengths. This study gives an insight of the 3D effects, i.e., beam-filling and horizontal photon transport effects, which contaminate cloud ice retrievals from both current and forthcoming satellite missions, with a focus on the Ice Cloud Imager (ICI). Accordingly, forward radiative transfer simulations were carried out for two incident angles at the ground, i.e., 53 • (slant view) and 0 • (nadir view); at two frequencies (186.3 GHz and 668 GHz); and with two footprint sizes, i.e., 6 km and 15 km, over several tropical and mid-latitude cloud scenes. The 3D synthetic scenes have been generated from two-dimensional (2D) CloudSat (Cloud Satellite) observations utilizing a stochastic approach.
By conducting a comprehensive statistical analysis over randomly sampled scenes, we demonstrated that, for mm/sub-mm wavelengths: 1. The horizontal photon transport effect is rather small, and thus, the errors induced by neglecting the 3D heterogeneities of clouds can be mostly linked to the beam-filling effect. In fact, conducting radiative transfer simulations by means of an independent beam approximation (IBA) induces a slight overestimation and chiefly random errors. Accordingly, the realistic but computational expensive 3D radiative transfer simulations could be excluded from retrievals of current and proposed satellite missions, providing the fact that the bias of the forward model is adjusted by a correction factor between 0.1 and 1 K. 2. On the contrary, performing one-dimensional (1D) radiative transfer simulations following the plane-parallel assumption divulges a significant beam-filling effect that increases primarily with frequency, and secondly, with footprint size and slant path. On this basis, we conclude that ice cloud retrievals that supply the optimal estimation approach (also known as 1D variational inversion algorithm, 1DVAR) by current or future satellite instruments are and will be adversely contaminated by the beam-filling effect. Consequently, this promotes the use of alternative methods, which are based on retrieval databases and utilize IBA simulations [92]. In fact, such a retrieval database has been suggested for the upcoming ICI operational product [10]. 3. In case IBA simulations are still not computationally feasible, then, among the two different 1D calculations modes, i.e., (a) HND-avg that employs the mean of the 3D hydrometeor number density field and (b) HC-avg that utilizes the 1D hydrometeor number density derived from the mean hydrometeor content, the latter one is preferable, since it leads to slightly smaller systematic errors-up to ≈1.9 K and a root mean square error (RMSE) of RMSE = 2.3 K at a frequency of 668 GHz with a footprint size of 15 km at which the largest errors were reported. 4. Our findings indicate that 1DVAR retrievals and data assimilation at sub-mm wavelengths will be largely affected by 3D effects. In our framework, a simple statistical correction scheme by means of a multiplication factor has been developed that compels the errors induced by neglecting domain heterogeneities to be more symmetric. This multiplication factor is the slope of the excellent linear correlation between 3D and 1D simulations. Under these circumstances, the RMSE of the 3D effect is reduced by up to 3.2 K. Note here that this was a first attempt towards a correction scheme for 1DVAR retrievals, with a potential use in all-sky assimilation; diminution of the 3D effect along the field of view of ICI at mm/sub-mm wavelengths instructs a separate study including all the channels of ICI, and thus, it was beyond the scope of this paper. 5. The considerable 3D effects reported here at sub-mm wavelengths with a footprint size of 15 km should be acknowledged for any future sensor design. Even though a radiometer operating at these wavelengths and characterized by smaller footprint size does not exist still, its cloud observations will be available in a finer spatial resolution and with a better overall precision considering the much smaller modeling uncertainties.
Future research could aim at further exploring the use of the current empirical correction scheme, or the development of a new empirical correction scheme for the beam-filling effect at mm/sub-mm wavelengths. In this work, due to the limitation in scattering databases that cover particle orientation at mm/sub-mm wavelengths, only randomly oriented hydrometeors have been considered. A subsequent step should be to include particle orientation in order to explore polarization effects that were introduced first by [44], and the recently published work by Brath et al. [93] makes such extension feasible.
Author Contributions: Both V.B. and P.E. conceived and refined the overall structure of the investigation. V.B. carried out the simulations, refined the data analysis, and wrote the manuscript. Both authors contributed to the interpretation of the results and to the manuscript's improvement. All authors have read and agreed to the published version of the manuscript.