Application of M5 Model Tree in Passive Remote Sensing of Thin Ice Cloud Microphysical Properties in Terahertz Region

: This article presents a new method for retrieving the Ice Water Path (IWP), the median volume equivalent sphere diameter (D me ) of thin ice clouds (IWP < 100 g/m 2 , D me < 80 µ m) in the Terahertz band. The upwelling brightness temperature depressions caused by the ice clouds at 325.15, 448.0, 664.0 and 874.0 GHz channels are simulated by the Atmospheric Radiative Transfer Simulator (ARTS). The simulated forward radiative transfer models are taken as historical data for the M5 model tree algorithm to construct a set of piecewise functions which represent the relation of simulated brightness temperature depressions and IWP. The inversion results are optimized by an empirical relation of the IWP and the D me for thin ice clouds which is summarized by previous studies. We inverse IWP and D me with the simulated brightness temperature and analyze the inversion performance of selected channels. The 874.4 ± 6.0 GHz channel provides the most accurate results, because of the strong brightness temperature response to the change of IWP in the forward radiative transfer model. In order to improve the thin ice clouds IWP and D me retrieval accuracy at the middle-high frequency channels in Terahertz band, a dual-channel inversion method was proposed that combines the 448.0 ± 3.0 GHz and 664.0 ± 4.2 GHz channel. The error analysis shows that the results of the 874.4 ± 6.0 GHz channel and the dual-channel inversion are reliable, and the IWP inversion results meet the error requirement range proposed by previous studies.


Introduction
Several researchers have studied the distribution characteristics and radiative impacts of thin ice clouds [1][2][3]. These researchers demonstrate that the existence of thin ice clouds is ubiquitous, and it plays an important role in future upper troposphere analysis. IWP, ice crystal shape and particle size are the critical parameters that describe and modify the net radiative effect of ice cloud [4]. However, the representation of the ice cloud in the General Circulation Models (GCMs) is currently incomplete [5,6], which may lead to producing errors in weather and climate forecasts. Improving the microphysical-property inversion accuracy of the thin ice clouds is conducive to studying the radiation effects of ice clouds and reducing the uncertainty in GCMs attributed to ice clouds. The Terahertz wave (0.1~10 THz, 30 µm~3 mm) detection is expected to bear a high potential to obtain ice clouds' microphysical properties [7]. The wavelength of THz radiation is in the middle of visible and microwave bands, comparable to the size of typical particles in ice clouds [8]. Therefore, THz radiation can well penetrate the cloud layer, and the changes of brightness temperature from ice clouds observed by the radiometers are correlated to ice mass.
THz wave radiometers and corresponding inversion algorithms have been developed to obtain the information of ice clouds [9]. IceCube and the Tropospheric Water and ICE CubeSat (TWICE) are the two small satellites launched by National Aeronautics and Space Administration (NASA) to study the ice cloud properties [10,11]. TWICE adopted the Bayesian Monte Carlo Integration (BMCI) method to retrieve the parameters of ice clouds [12]. The Bayesian methodology represents the information from prior knowledge, observations, and the forward model probabilistically. If the prior information and the physical assumptions are insufficient, the inversion could be unsuccessful for some measurements [13]. A neural networks (NNs) method is used in the International Submillimeter Airborne Radiometer (ISMAR) to retrieve the quantities of ice clouds [14]. ISMAR is funded by UK Meteorological Office and European Space Agency (ESA), serves as an airborne demonstrator for the upcoming Ice Cloud Imager (ICI) [15]. The NNs inversion method takes the neuron as the basic processing unit to set up a nonlinear relationship between brightness temperature depression and parameters of ice clouds with a training dataset simulated by ARTS [16]. All the assumptions on the retrieval method are condensed in the database for NNs and the BMCI method, and the accuracy of the method has relied on the quality of the inversion priori database. A lookup table, optimal estimation and the multichannel linear regression method are also used to retrieve the parameters of ice clouds. Applications of optimal estimation method are infrequent, mainly because the assumptions of this method ignore the correlations between variables and the loop iterations in the algorithm, and reduce the computational efficiency [17]. The two-dimensional multiple lookup table method studies the relationships between the parameters such as brightness temperature depression, bright temperature depression slope, IWP and D me [18]. The method optimizes the relations by interpolation and normalization to set up the lookup table at typical channels. Multichannel linear regression is a relatively efficient method to retrieve IWP [19], but this method and the lookup table method have strict requirements of the channel combinations which may limit their applications.
This work intends to develop a less time-consuming but accurate inversion method for passive remote sensing ice clouds in the THz band. We apply the M5 algorithm to retrieve the IWP and the D me of thin ice clouds. The algorithm sets up a multivariate linear model to fit the nonlinear relationship between IWP and brightness temperature depression, and adopts an empirical relation to restrain and optimize the inversion results. This less time-consuming method can be applied in a single-channel inversion to obtain both the IWP and the D me of ice clouds. The parameters of the forward radiative transfer model and the detail inversion process are described in Section 2. In Section 3, we discuss the inversion results from single -channel and the dual-channel algorithms with simulated brightness temperature depression. The inversion method summaries and inversion error analysis are presented in the last section.

Forward Radiative Transfer Model
The Atmospheric Radiative Transfer Simulator (https://www.radiativetransfer.org/ accessed on 25 May 2021) is a public domain software for radiative transfer simulations in the thermal spectral range (microwave to infrared) [20]. Several intercomparison studies have validated the models calculated by ARTS [21,22]. In this work, we simulate the brightness temperature depression (∆T) of upwelling radiation caused by the ice clouds at 325.15, 448.0, 664 and 874.4 GHz. These channels are selected at the atmosphere windows or gas absorption regions suggested by Liu et al [23]. The inputs of the ARTS include atmospheric profiles, particle absorbing and scattering attenuation data. Atmospheric profiles such as temperature, relative humidity, pressure and trace gaseous composition are derived from the Fast Atmospheric Signature CODE (FASCOD) standard atmospheric profile database [24]. The absorbing attenuation spectral parameters of gas molecules in the atmosphere are extracted from the HITRAN2012 database [25]. Gang et al. [26] calculated a single scattering property dataset of non-spherical particles at 90-874 GHz with the discrete-dipole approximation (DDA) algorithm. The scattering attenuation parameters of column particles are derived from this library, and the Discrete Coordinate Iterative (DOIT) method of ARTS is used to simulate the scattering process of cloud particles.
The characteristics of ice clouds microphysical properties are summarized from in situ measurements with aircraft-borne sensors [27]. The effects of these properties on THz radiation have also been investigated in previous studies [19,28]. Normally, the maximum diameter of a thin ice cloud particle is between 30 and 140 µm. We convert the maximum diameter into D me by the fitted relations which conclude from a large number of field experiments [29]. Accordingly, the brightness temperature depressions caused by the ice clouds with D me in 10-80 µm, IWP in 0-1000 g/m 2 are simulated. Figure 1 presents the simulation of a homogeneous ice cloud at an altitude of 10.5 km in the midlatitude winter using the ARTS forward radiative transfer model. The figure shows that the 874 GHz channel are sensitive to the change of IWP, the max brightness temperature depression (the max ∆T) at 874.4 ± 6.0 GHz channel is above 50 K while the value is under 20 K for 325.15 ± 9.5 GHz channel.

Building an M5 Model Tree
The M5 scheme and construction have been described by Quinlan et al. [30], which can be used to predict the continuous response variable values. It is a binary regression tree with a linear regression function at terminal (leaf) nodes. Classification and Regression Trees (CART) use historical data to construct a decision tree to classify the new data [31]. It arranges a set of tests to find the best split to divide the data into two segments with maximum homogeneity. Model trees have a similar structure to the regression trees, but they use the fitted linear regression model to predict the response values. Several studies have applied the M5 model tree in the fields of atmosphere science and hydrology [32,33], and the comparison between the model tree and neural network method for predictions. The results from the comparison showed that the model trees were marginally more accurate than NNs and less time-consuming [32].
There are two stages to generate a model tree. The first stage includes splitting the instance space recursively to construct a regression tree, in which the maximization of the intra-subset variation of the target value is used as the splitting criterion. The variation is measured by the standard deviation of the values that reach the node from the root to the branch with calculating the expected reduction in errors as a result of testing each attribute at the node [30], and the split which maximizes the expected error reduction is chosen. To prevent overfitting, splitting stops if only a few instances remain. The length of the instance is termed as a pre-pruning condition. After the full model tree is constructed, the tree is pruned by replacing subtrees with linear regression functions if they lower the estimated error. The test data are randomly selected from the historical data, and then the errors of the values predicted by the full model tree and the pruned model tree, respectively. If the pruned model tree decreases the estimated error, the subtree is replaced. We repeat the prune process and implement 10 cross-validations to construct the final model tree. The full process of building a model tree is shown in Figure 2. In this study, the whole forward radiative models are used to construct model trees for each selected channel. The feature of the historical data is IWP and the response variable is the brightness temperature depression caused by the simulated ice clouds, and the D me is noted as a label for each fitted piecewise function.

Restrain Inversion Result with an Empirical Relation
General circulation models predict D me from relationships between the cloud temperature and Ice Water Content (IWC) [34]. A previous study verified that a global parameterization may be suitable for GCMs and parameterizes D me as functions of IWC and IWP by analyzing the collocated Radar-Lidar Geometrical Profiling (GEOPROF) data [35]. The relationship of IWP and D me is described by the expression: where α 0 = 9.945, α 1 = −11.245, α 2 = 27.426, α 3 = −9.192 and α 4 = 0.974. By assuming a vertically constant IWC which is in accordance with the calculated forward radiative transfer model, IWC is determined as a ratio between IWP and the vertical extent of the cloud layer. Thus, D me can be parameterized as a function of the IWP. The relationship between IWP and D me is used to optimize the retrieved results. As described in Section 2.2, after fitting the forward radiation transfer model with the model tree algorithm, a set of piecewise functions of IWP and ∆T with labeled D me are generated. Partial fitted results at 664.0 ± 4.2 GHz are presented in Table 1 as an example. The piecewise function of each labeled D me is listed in the corresponding column. IWP can be solved for each labeled D me with Formula 1. To compensate for the errors in the detection and inversion process, IWP in the range of IWP ± IWP × 0.15 are accepted in the inversion. The Formula 1 specifies acceptable value range of IWP for the sequence of D me . The subsection in the piecewise function that covers the specified IWP range in the forward radiative model is noted for inversion. Each noted equation related to a determined D me , representing a fitted relationship of simulated IWP and ∆T. The selected equations of the piecewise functions at 664.0 ± 4.2 GHz are shown in Table 2. To retrieve the ice cloud microphysical properties with the detected brightness temperature depressions, all the selected functions of labeled D me are used to retrieve IWP, then a set of IWP and corresponding D me can be determined. The retrieved IWP-D me pair which fits the empirical relationship in Formula 1 is selected. If the inversed IWP is within the acceptable range of IWP' associated with D me , then the IWP and corresponding D me are the optimal inversion results. Figure 3 illustrates the retrieving process with the fitted relationship at 664.0 GHz. The grey rectangle near the x-axis is denoting the accepted IWP range for simulated ice clouds which correspond to D me = 73 µm. The three solid black lines represent the selected equations of 73, 65 and 54 µm that cover the conformable part of the forward radiative transfer model (The conformable part is represented by the solid red line).
As shown in Figure 3, the ∆T for illustrating the retrieving process is 15K, the IWP are retrieved with the equations that labeled D me is 73, 65 and 54 µm respectively. There is one retrieved IWP in the IWP' range which is specified by 73 µm D me , while the IWP retrieved by the other two equations deviate from the corresponding accepted range. In this way the optimal inversion result is determined. The flow chart of inversion process is presented in Figure 4.

Application of Inversion Method and Discussion
The retrieval performance of the method is discussed in this section. Linear interpolation is performed for the forward radiative transfer model to increase the

Application of Inversion Method and Discussion
The retrieval performance of the method is discussed in this section. Linear interpolation is performed for the forward radiative transfer model to increase the resolutions of IWP and D me before inversions. To evaluate the inversion results, part of the forward radiative transfer model where IWP and D me conform to Formula 1 is noted as testing data, and the IWP and D me in the test set are taken as the true value for inversion. After analyzing the instrument noise levels in the existing Terahertz radiometers [15,36], a random noise in the range of ±1 K is introduced to ∆T in the test data. We perform 10 times inversion for each simulated case to improve the reliability of the results.

Results from Single-Channel Inversion
The inversed IWP and D me of the 325.15, 448.0, 664.0 and 874.0 GHz single-channel inversion are presented in Figures 5 and 6, respectively. The dotted lines in Figures 5 and 6 indicate a ±20% relative accuracy range for inversion results. In both figures, the inversed results show an increasing consistency between the simulated true values and the inversed values with an increase of the channel frequency. While increasing the frequency of the selected channels, the response ∆T in the forward radiative transfer model to the change of IWP becomes stronger. This feature seems to implicate that increasing the frequency of the channel can improve the inversion accuracy. Overall, 874.4 ± 6.0 GHz is the most ideal single channel for this method to retrieve IWP and D me of thin ice clouds. Inversion absolute error of 325 GHz channel is 27.2 g/m 2 for IWP, 14.28 µm for D me , the inversion results are affected by the introduced noise. A primary reason for the unappealing results is that the ∆T is insensitive to IWP at 325 GHz channel. The inversion results for low IWP (<40 g/m 2 ) from 448.0 GHz channel are also inaccurate, but the relative inversion errors are mostly within 20% for the retrieved IWP that is above 50 g/m 2 . The 664.0 and the 874.0 GHz channels can provide relatively accurate inversion results to other selected channels, especially for the low IWP. However, after performing interpolation for the forward radiative model, the acceptable range for IWP might overlap for several labeled D me , which may produce errors in the D me inversion process. It can be seen from Figure 6c,d that the errors are more obvious for high-frequency channels for retrieving high IWP (>60 g/m 2 ), mainly because the inversion deviations of labeled D me would bring more errors to IWP inversion for channels with higher frequency in THz band. As a consequence, there are some deviations for high-value IWP retrieved by 664.0 and 874.0 GHz single-channel algorithm.

Results from Dual-Channel Inversion
To obtain the accurate inversion results from the middle-high frequency channels in the THz band, a dual-channel inversion scheme is proposed that combines 448.0 ± 3.0 GHz and 664.0 ± 4.2 GHz channels to retrieve IWP and D me jointly.
In the inversion process, IWP is firstly retrieved by these two channels. If the IWP inversed by 664.0 ± 4.2 GHz is below 40 g/m 2 , the inversed IWP and D me are taken as final inversion results. Whereas, if the inversed IWP is above 40 g/m 2 , we first examine D me inversion results from these two channels. If retrieved D me are identical, then the final IWP inversion results are determined by calculating the mean value of IWP retrieved by both channels; otherwise, the inversion uses D me retrieved from 664.0 ± 4.2 GHz but the labeled relations fitted from 448.0 ± 3.0 GHz to inverse IWP as a final result. This way reduces the errors caused by retrieved D me deviation of 664.0 ± 4.2 GHz. Figure 7 presents the retrieved IWP, D me from a dual-channel inversion. The mean inversion absolute errors of IWP and D me are decreased; the mean absolute errors reduce to 7.930 g/m 2 for IWP and 5.997 µm for D me . The results reveal that the utilization of a dual-channel inversion improves the performance of 448.0 ± 3.0 or 664.0 ± 4.2 GHz channel inversion as well.

Results from Dual-Channel Inversion
To obtain the accurate inversion results from the middle-high frequency chann the THz band, a dual-channel inversion scheme is proposed that combines 448.0 ± 3.0 and 664.0 ± 4.2 GHz channels to retrieve IWP and Dme jointly.
In the inversion process, IWP is firstly retrieved by these two channels. If the inversed by 664.0 ± 4.2 GHz is below 40 g/m 2 , the inversed IWP and Dme are taken as inversion results. Whereas, if the inversed IWP is above 40 g/m 2 , we first examine inversion results from these two channels. If retrieved Dme are identical, then the final inversion results are determined by calculating the mean value of IWP retrieved by channels; otherwise, the inversion uses Dme retrieved from 664.0 ± 4.2 GHz but the lab relations fitted from 448.

IWP Inversion Error Analysis
This section presents further examination of 664.0 ± 4.2 and 874.4 ± 6.0 GHz singlechannel and the dual-channel IWP inversion results. We summarize the optimal inversion channel for the method and the expected errors of each selected inversion channel.
Previous studies proposed a requirement range for an absolute error and a relative error to assess the accuracy of the inversions [37]. Neither absolute nor relative errors alone are adequate to describe a real measurement system on account of the large dynamic range of the inversed quantities. The absolute error requirement range for IWP below 20 is 1-10 g/m 2 , for IWP which are above 20 g/m 2 , the relative error requirement range is 10-50%. Accordingly, the absolute errors are calculated for inversed IWP that is below 20 g/m 2 , otherwise the relative errors are applied to assess the inversion accuracy. The inversed IWP absolute error and relative error of 664, 874 GHz and the dual-channel algorithm are displayed in Figure 8a,b. The boxplots show the quartile, median, maxima and minima of the error information, and the data that passed the low and high extremum are defined as abnormal values. The upper and lower extremum are calculated by: where qua up and qua lo are the upper and the lower quartile of error data, respectively.

Conclusions
This article describes a new method for retrieving IWP and D me of the thin ice clouds in Terahertz region. The method applies the M5 model tree algorithm in the singlechannel and dual-channel inversion. We perform inversions with brightness temperature simulated by ARTS for the frequency of 325.15, 448.0, 646.0, 874.4 GHz. Then, the inversion performance at each channel is summarized. In the inversion process, a random noise within 1 K was added to the simulated brightness temperature depression to verify the effectiveness of the method.
The single-channel inversion results error analysis demonstrates that the 874.4 GHz channel is an optimal inversion channel for this method. The inversion results of the 325.15 and 448.0 GHz channels are not ideal. However, the dual-channel inversion manages to decrease inversion errors for the middle-high frequency channels by combining 448.0 and 646.0 GHz channels. The overall error analysis shows that the inversion accuracy of 874.4 GHz single-channel and the dual-channel inversions meet the requirement range proposed by previous studies. The essential idea of this article is based on applying the model tree algorithm in constructing piecewise functions for forward radiative transfer models, and restraining the inversion results with an empirical relation between the predictor variables. One of the merits of the model tree algorithm is that it can process the high dimensional data efficiently. This study uses the model trees to construct linear equations with IWP and brightness temperature. To fully exploit the advantage of the method, more features of the ice clouds such as the cloud height and the vertical extent of cloud could be taken into consideration in the future works. Additionally, the accuracy of the empirical relation used in the inversion has an evident impact on the retrieved results. Thus, to summarize overall and practical empirical relations for inversion is another important work.