Rainfall Estimation from Tropospheric Attenuation Affecting Satellite Links

A novel methodology for estimating rainfall rate from satellite signals is presented. The proposed inversion algorithm yields rain rate estimates by making opportunistic use of the downlink signal and exploiting local ancillary meteorological information (0 ◦C isotherm height and monthly convectivity index), which can be extracted on a Global basis from Numerical Weather Prediction (NWP) products. The methodology includes different expressions to take the different impact of stratiform and convective rain events on the link into due account. The model accuracy in predicting the rain rate is assessed (and compared to the one of other models), both on a statistical and on an instantaneous basis, by exploiting a full year of data collected in Milan, in the framework of the Alphasat Aldo Paraboni propagation experiment.


Introduction
The observation, monitoring, and estimation of rainfall are matters of major concern for a number of reasons. For instance, public administrators aim at guaranteeing an adequate level of security for people living and operating on a given territory: for this purpose, rainfall forecasting is key for the assessment of hydrogeological risks due to landslides and floods. An accurate measurement of total rainfall amount is also critical in climatological studies, since it represents one of the main climate indices. Moreover, the knowledge of the real-time rain rate during precipitation events is useful in urban areas to estimate water runoff from the drainage system.
A very high spatial and temporal variability characterizes rainfall [1][2][3]: it is then difficult to obtain accurate measurements with proper spatial density and temporal detail. Several instruments and techniques have been developed to gauge this quantity (e.g., the use of weather radars, rain gauges, infrared satellite imagery, microwave links, etc.). Each of them presents advantages and limitations [4]: for example, remote sensing measurements aim at estimating the precipitation field over wide areas, but typically with limited spatial resolution and accuracy; on the contrary, point observations, e.g., collected by rain gauges and disdrometers, deliver precise local measurements, but they hardly provide large and dense coverage.
Recently, increasing attention has been paid to the exploitation of satellite downlink signals for the estimation of rainfall, as most of the commercial satellite communication systems currently operate in the Ku band (12)(13)(14)(15)(16)(17)(18) and the K band (18)(19)(20)(21)(22)(23)(24)(25)(26): in fact, in such a frequency range, rain significantly impairs the propagation of electromagnetic waves [5,6]. Some approaches have been developed for estimating the intensity of precipitation on a large scale from space-borne signals, with a spatiotemporal resolution between that of rain gauges and of weather radars. In fact, the opportunistic Information 2020, 11, 11 2 of 19 use of satellite downlink signals to estimate rainfall combines the benefit of covering areas where ground instrumentation or radar networks are not available, and the requirement of high temporal resolution measurements, which is in the order of one minute [7]: this approach would allow providing two-dimensional maps of the rain intensity, even across continental regions (e.g., Europe) by using a sufficiently dense network of ground receivers.
The methodologies that are currently available to assess rainfall rate via satellite links are of limited accuracy and reliability and, in some cases, of complex implementation [8][9][10][11]: as a consequence, it becomes necessary to investigate alternative estimation methods, which only require few measurable physical parameters as input.
This contribution presents a novel estimation methodology that aims to estimate the rain rate by opportunistically exploiting local downlink satellite communication signals and by making use of local ancillary meteorological data, namely the hourly zero-degree isotherm height and the monthly convectivity index. The model also includes different expressions to take the different impact of stratiform and convective rain events on the link into due account. More specifically, the remainder of the paper is structured as follows. Section 2 presents the Earth-space propagation experiment carried out at Politecnico di Milano, along with the experimental setup used in this work. Section 3 describes the experimental datasets. Section 4 describes the proposed prediction methodology in detail, while Section 5 summarizes the full procedure for rain rate estimation. Lastly, Section 6 discusses the performance of the proposed model (as well as of other methodologies that are available in the literature) and Section 7 draws some conclusions.

Experimental Setup
A relevant part of the data used during this work is collected within the measurement campaign of the Alphasat Aldo Paraboni propagation experiment [12], whose main goal is to better investigate the propagation impairments that affect Ka/Q/V band satellite channels. In 2014, NASA Glenn Research Center (GRC) joined the experimental campaign by installing the equipment to receive the beacon signals transmitted by the Aldo Paraboni Payload (also known as Technology Development Payload 5-TDP 5) [13] aboard the Alphasat geosynchronous Earth orbit satellite at Politecnico di Milano. The payload is equipped with two beacons that transmit coherent, continuous-wave signals at 19.701 and 39.402 GHz, respectively.
The NASA experimental equipment is installed on the rooftop of the Dipartimento di Elettronica, Informazione e Bioingegneria (DEIB) at Politecnico di Milano, Italy. The instruments, as depicted in Figure 1, include:

Experimental Datasets
The model developed in this work makes use of the following data: • Ka-band rain attenuation time series derived from the received beacon power levels.

Optical Weather
Ka-band Q-band  Two beacon receivers (Ka band and Q band), each of which consists in a Cassegrain antenna (diameter of 0.6 m and 1.2 m for the Ka and Q band, respectively), and a down-conversion stage enclosed in a box positioned just behind the antenna. The received power level is sampled with a rate of 8 Hz.

•
A Thies Clima Laser Precipitation Monitor, an optical disdrometer counting and classifying hydrometeors according to their size and falling velocity, from which the rain rate R can be calculated with 1-min integration time.

•
A weather station, monitoring wind speed and direction, pressure, temperature, and relative humidity, with a sampling interval of 1-min.

•
A tipping bucket rain gauge, providing an additional measurement of R (1-min integration time), though with a more limited accuracy if compared to the disdrometer.

Experimental Datasets
The model developed in this work makes use of the following data: • Ka-band rain attenuation time series derived from the received beacon power levels.

•
Rain rate time series and specific attenuation of the site calculated from the drop size distribution (DSD) that was measured by the disdrometer.

•
The zero-degree isotherm height, h 0 , extracted from the vertical temperature profiles produced by the European Centre for Medium-Range Weather Forecasts (ECMWF).
The measurement period spans from 1st January 2017 to 31st December 2017. A preliminary selection was performed to rule out rain events when precipitation occurred along the link, but not directly over the receiver site, as well as days when at least one of the two instruments (disdrometer and the beacon receiver) was not operational. Overall, the rain events that were used in this work were collected in 70 different days.

Rain Attenuation Time Series
Rain-induced attenuation A R is extracted starting from the received beacon power levels that were collected by the NASA Ka-band receiver located in Milan. The typical, well-established approach to derive A R from the received beacon power P R is to first low-pass filter P R with a typical cut-off frequency of 0.03 Hz, to remove scintillations [14]. Afterwards, rain events are identified, usually both by taking advantage of the local rain sensors (if present) and by inspecting the trend of P R : indeed, the impact of rain on the link might be longer than what is recorded by the disdrometer, especially at the beginning and/or the end of the event. Afterwards, rain attenuation is calculated by subtracting from P R the power level that is the linear interpolation of P R just before the beginning and after the end of each event [14]. Lastly, the attenuation data are averaged to obtain a sampling period equal to the one of the disdrometer (1-min per sample).
An example of such a well-established procedure, which, for instance, is also employed in [14], is depicted in Figures 2-4, where the sequence of operations (low-pass filtering, rain event detection, and estimation of the reference power level during the rain event through linear interpolation) are applied to the received power on 7th September 2017. Figure 5 plots the resulting rain attenuation time series.

Rain Rate Time Series and Specific Attenuation of the Site
A LabVIEW software was designed at NASA GRC to monitor and process the data that were obtained by the Thies Clima Laser Precipitation Monitor. The measured rain rate time series and the DSD data are retrieved from the output files provided by such a software. In particular, the drop size distribution ( ) [15], which represents the number of particles per cubic meter having diameter , can be used to derive the specific attenuation of the site, [16]: where is the frequency of interest and is the extinction section of the raindrops with diameter D, which were computed while using the T-Matrix method [17].
The regression coefficients and (reported in Table 1) were derived by fitting the power law relationship γ = kR α [16] on the points of the scatter plot shown in Figure 6, where the 1-min values of the specific rain attenuation are plotted against the rain rate with the same integration time. Figure 6 also shows the specific rain attenuation calculated while using recommendation ITU-R P.838-3 [18] (the corresponding and coefficients are reported in Table 1), which tends to overestimate the specific attenuation derived from the disdrometer in Milan.

Rain Rate Time Series and Specific Attenuation of the Site
A LabVIEW software was designed at NASA GRC to monitor and process the data that were obtained by the Thies Clima Laser Precipitation Monitor. The measured rain rate time series and the DSD data are retrieved from the output files provided by such a software. In particular, the drop size distribution N(D) [15], which represents the number of particles per cubic meter having diameter D, can be used to derive the specific attenuation of the site, γ R [16]: where f is the frequency of interest and C ext is the extinction section of the raindrops with diameter D, which were computed while using the T-Matrix method [17]. The regression coefficients k and α (reported in Table 1) were derived by fitting the power law relationship γ = kR α [16] on the points of the scatter plot shown in Figure 6, where the 1-min values of the specific rain attenuation are plotted against the rain rate with the same integration time. Figure 6 also shows the specific rain attenuation calculated while using recommendation ITU-R P.838-3 [18] (the corresponding k and α coefficients are reported in Table 1), which tends to overestimate the specific attenuation derived from the disdrometer in Milan.

Rain Rate Time Series and Specific Attenuation of the Site
A LabVIEW software was designed at NASA GRC to monitor and process the data that were obtained by the Thies Clima Laser Precipitation Monitor. The measured rain rate time series and the DSD data are retrieved from the output files provided by such a software. In particular, the drop size distribution ( ) [15], which represents the number of particles per cubic meter having diameter , can be used to derive the specific attenuation of the site, [16]: where is the frequency of interest and is the extinction section of the raindrops with diameter D, which were computed while using the T-Matrix method [17].
The regression coefficients and (reported in Table 1) were derived by fitting the power law relationship γ = kR α [16] on the points of the scatter plot shown in Figure 6, where the 1-min values of the specific rain attenuation are plotted against the rain rate with the same integration time. Figure 6 also shows the specific rain attenuation calculated while using recommendation ITU-R P.838-3 [18] (the corresponding and coefficients are reported in Table 1), which tends to overestimate the specific attenuation derived from the disdrometer in Milan.

Zero-Degree Isotherm Height
The freezing level corresponds to the zero-degree isotherm height. Below such altitude, the temperature is above the freezing point and rainfall may occur. The knowledge of the freezing level is necessary for the estimation of the equivalent rain height h R , in turn used in any rain attenuation prediction model: the latter quantity represents the altitude of an equivalent rain layer, which is typically higher than the freezing level and it depends on the different precipitation type (stratiform and convective).
The zero-degree isotherm height time series h 0 (t) for Milan are derived from ECMWF high-resolution vertical profiles of temperature. Specifically, the data are extracted from the ERA5 database, whose temporal resolution is one hour (from 00:00 to 23:00 UTC), with a spatial resolution of 0.25 • × 0.25 • in latitude and longitude. At each time instant t, the temperature profile provides the temperature T(h) as a function of the altitude h; h 0 (t) is then derived by linearly interpolating the two h values closer to T(h) = 0 • C.

The Rainfall Rate Estimation Model
The proposed rain rate estimation algorithm is a typical inverse model; specifically, it is a closed-form relationship that expresses the rain rate as a function of two physical quantities, namely the rain attenuation A R (obtained from the link) and the equivalent rain height h R (estimated from ECMWF data, as described more in detail in Section 4.3).
The chosen model to be inverted is the prediction method that is described in Section 2.2.1.1 of recommendation ITU-R P.618 [19], which provides as output the rain attenuation A R,p exceeded for an arbitrary time percentage p in a year, as a function g(R, h R , k, α) of the rain rate R (mm/h), of the equivalent rain height h R (km), and of the k and α power-law coefficients of Table 1: (2) The model, whose full mathematical framework can be found in [19] (and is not reported here for the sake of brevity), has proven to reliably predict rain attenuation on a statistical basis [20][21][22]: the idea that is put forth in this contribution (and verified in the remainder of the contribution) is that such a model can also be used to provide accurate results on an instantaneous basis. It is worth underlying that part of the good prediction accuracy of this model is ascribable to the inclusion of path reduction factors, analytical expressions that allow for taking the inhomogeneity of precipitation along the link in due account. As a result, the model offers an accurate prediction of the rain attenuation along the path from the sole knowledge of the rain rate measured at the receiver site; this means that it is possible to predict point rain rate values to be compared against reference precipitation data collected at the receiver site through the inversion of the model.
As is clear from [19], the function g includes various non-linear relationships between A R and R: as a consequence, no analytical inversion for g(R, h R , k, α) exists.

Approximate Solution to the Inversion Problem
When considering that an analytical inversion of g is not possible, the rain rate can be alternatively expressed as: where q is an approximation of g −1 .
Information 2020, 11, 11 7 of 19 After investigating the trend of (2) conditioned to a given value of h R , and while using the specific power law coefficients of the site (listed in the third column of Table 1), it turns out that g −1 can be accurately approximated by the following third order polynomial q, which is a function of R and h R : The three regression coefficients a, b and c are the following functions of h R : Figure 7 compares g −1 , as calculated by implementing the expressions in [19], and q, obtained by using Equations (4)-(7), for different values of h R . Table 2 lists the indicators that were used to evaluate the fitting accuracy of q. In particular, we have calculated the residual ε =x − x (x and x represent the values of the fitting curve and of the input data points, respectively), and the following statistical indicators: the average (E), the R-Square (R 2 ), and the Root Mean Square (RMS). The fitness accuracy values in Table 2 were obtained by averaging results over 200 different curves (as an example, Figure 7 shows four of them), to consider a broad range of possible h R values, from 1.65 km up to 8 km. Regarding the coefficients a, b, and c, as the E and RMS values of the fitting errors are infinitesimal, they have been rounded to zero. can be accurately approximated by the following third order polynomial , which is a function of and ℎ : The three regression coefficients , and are the following functions of ℎ : (ℎ ) = 5.7255 ⋅ ℎ . + 0.9172. (7) Figure 7 compares g −1 , as calculated by implementing the expressions in [19], and q, obtained by using Equations (4)-(7), for different values of ℎ . Table 2 lists the indicators that were used to evaluate the fitting accuracy of q. In particular, we have calculated the residual = − ( and x represent the values of the fitting curve and of the input data points, respectively), and the following statistical indicators: the average (E), the R-Square ( ), and the Root Mean Square (RMS). The fitness accuracy values in Table 2 were obtained by averaging results over 200 different curves (as an example, Figure 7 shows four of them), to consider a broad range of possible ℎ values, from 1.65 km up to 8 km. Regarding the coefficients a, b, and c, as the E and RMS values of the fitting errors are infinitesimal, they have been rounded to zero.  Table 1.

Equivalent Rain Height Model
Stratiform and convective rain events have a different impact on satellite links, not only in terms of specific attenuation, but mainly as for the equivalent rain height, which can change on an event basis. In particular, the proposed model includes two distinct formulations for ℎ :  Table 1.

Equivalent Rain Height Model
Stratiform and convective rain events have a different impact on satellite links, not only in terms of specific attenuation, but mainly as for the equivalent rain height, which can change on an event basis. In particular, the proposed model includes two distinct formulations for h R : • In the case of stratiform events, the presence of the melting layer, which is also known as bright band, induces extra attenuation on electromagnetic waves propagating through the troposphere: as a matter of fact, in this case, the rain rate is constant approximately up to the freezing level, below which frozen precipitation particles start to melt during their fall (melting layer). According to the SC EXCELL model [23], the additional attenuation due to the melting layer is taken into account by means of an equivalent rain height h BB , to be added to h R , with the latter being usually set equal to the zero-degree isotherm height h 0 . In mathematical terms, the stratiform rain height is given by [23]: where f is the frequency in GHz.

•
In the case of convective events, the corresponding rain height is higher than the zero-degree isotherm height because of the strong updrafts and downdrafts, which push water particles up beyond h 0 and also prevent the formation of a definite melting layer. Consequently, the convective rain height is an increased value of h 0 : In (10), the rain height enhancement factor τ(m) is higher than 1 and depends on m, namely the month number of the calendar year.

Extraction of the 0 • C Isotherm Height
As shown in (8) and (10), the equivalent rain height h R depends on the "instantaneous" value of the 0 • C isotherm height h 0 , which, as mentioned in Section 3.3, is extracted from the ERA5 database with a time resolution of an hour. Specifically, each sample of h 0 (t * ) is used in the model for the 60 min following t * .

Stratiform-Convective Discrimination Method
A key step of the model consists in classifying rain events as stratiform or as convective, hence in deciding whether (8) or (10) should be used to calculate the rain height. In temperate regions (like in Italy), there is a clear prevalence of stratiform events during cold months, and of convective ones during warm/hot months. This allows for adopting a monthly-based approach (rather than focusing on single events) to classify rain events and to, consequently, determine the values of τ(m). To this purpose, as τ(m) is expected to be tightly linked to the convectivity level of each month, we have investigated the relationship between τ(m) and β(m), with the latter being defined as the monthly ratio between the local convective and total accumulated rain amounts, R T,cnv (mm) and R T (mm), respectively: Figure 8 shows the yearly trend of β(m) for Milan, as derived from the ERA15 database that was produced by the ECMWF. The figure also reports the threshold on β(m), set to distinguishing between stratiform and convective months (black dashed line): we have chosen 0.1 based on the observation of an extensive database of rain events recorded in Milan. Using that threshold, six months of the year are classified as stratiform, while the remaining ones are labelled as convective. Therefore:

Derivation of the Rain Height Enhancement Factor
The values of ( ) were determined by maximizing the agreement between the annual Complementary Cumulative Distribution Function (CCDF) of the rain rate, P(R), as measured by the disdrometer and as estimated while using the rain rate prediction model; the two curves are reported in Figure 9 using a red and a blue line, respectively. We have obtained a first guess of the coefficients, named ̂ ( ) , which are listed in the first row of Table 3, as a result of this maximization process. It should be noted that the value of ̂ ( ) associated to August ( = 8) is missing, as no propagation data are available due to some issues with the beacon receiver.
We have investigated the relationship between ( ) and ( ) for the months classified as convective (that is for 4 ≤ ≤ 9) with the aim of devising a model with a broad applicability. Figure 10 shows that the two quantities follow the same trend, which confirms that the higher the convectivity level of a given month, the higher the increase in the rain height will be. In particular, the product ( ) between ( ) and ( ) turns out to be well approximated by the following third order polynomial, which is only function of : As a result, the rain height enhancement coefficient can be calculated starting from the local values of ( ), which can be obtained on a Global basis, as: The resulting values of ̂ ( ) are reported in the second row of Table 3. Figure 9 instead reported the annual P(R) estimated while using the coefficients ̂ ( ) (green line), and compared to the other two curves mentioned above: the agreement between the P(R)s estimated using ̂ ( ) and ̂ ( ) is very good. Table 4 also confirms this result, which lists the mean and RMS values of the estimation error ( ), in turn defined as the difference, for each probability level , between the measured (disdrometer) and estimated (model) rain rate.

Derivation of the Rain Height Enhancement Factor
The values of τ(m) were determined by maximizing the agreement between the annual Complementary Cumulative Distribution Function (CCDF) of the rain rate, P(R), as measured by the disdrometer and as estimated while using the rain rate prediction model; the two curves are reported in Figure 9 using a red and a blue line, respectively. We have obtained a first guess of the coefficients, namedτ 0 (m), which are listed in the first row of Table 3, as a result of this maximization process. It should be noted that the value ofτ 0 (m) associated to August (m = 8) is missing, as no propagation data are available due to some issues with the beacon receiver.
We have investigated the relationship between τ(m) and β(m) for the months classified as convective (that is for 4 ≤ m ≤ 9) with the aim of devising a model with a broad applicability. Figure 10 shows that the two quantities follow the same trend, which confirms that the higher the convectivity level of a given month, the higher the increase in the rain height will be. In particular, the product f p (m) between τ(m) and β(m) turns out to be well approximated by the following third order polynomial, which is only function of m: As a result, the rain height enhancement coefficient can be calculated starting from the local values of β(m), which can be obtained on a Global basis, as: The resulting values ofτ p (m) are reported in the second row of Table 3. Figure 9 instead reported the annual P(R) estimated while using the coefficientsτ p (m) (green line), and compared to the other two curves mentioned above: the agreement between the P(R)s estimated usingτ 0 (m) andτ p (m) is very good. Table 4 also confirms this result, which lists the mean and RMS values of the estimation error r (P), in turn defined as the difference, for each probability level P, between the measured (disdrometer) and estimated (model) rain rate. Information 2018, 9, x 10 of 18

Full Procedure for Rain Rate Estimation
The full procedure to estimate the rain rate consists in the following steps: (1) measure the instantaneous rain attenuation affecting the satellite link, ( ) [dB]; (2) obtain the concurrent value of the freezing level, ℎ ( ) [km]. This quantity can be extracted, for example, from vertical temperature profiles, as explained in Section 3.3; (3) if 4 ≤ ≤ 9: (a) compute ( ) from (11);

Full Procedure for Rain Rate Estimation
The full procedure to estimate the rain rate consists in the following steps: (1) measure the instantaneous rain attenuation affecting the satellite link, ( ) [dB]; (2) obtain the concurrent value of the freezing level, ℎ ( ) [km]. This quantity can be extracted, for example, from vertical temperature profiles, as explained in Section 3.3; (3) if 4 ≤ ≤ 9: (a) compute ( ) from (11);   Table 4. Mean and RMS values for the comparison between the curves reported in Figure 9.

Full Procedure for Rain Rate Estimation
The full procedure to estimate the rain rate consists in the following steps: (1) measure the instantaneous rain attenuation affecting the satellite link, A R (t) [dB]; (2) obtain the concurrent value of the freezing level, h 0 (t) [km]. This quantity can be extracted, for example, from vertical temperature profiles, as explained in Section 3.3; (3) if 4 ≤ m ≤ 9: (a) compute β(m) from (11); (b) evaluateτ p (m) using (14); (c) derive the equivalent rain height as (4) Calculate the coefficients a(h R ), b(h R ) and c(h R ), using (5), (6), and (7), respectively. (5) Estimate the rain rate as:

Performance Evaluation
The calibration of the model's parameters was carried out on a statistical basis, i.e., by maximizing the agreement between the measured and estimated P(R)s. However, it is worth verifying the accuracy in estimating the rain rate applicability when the model is used on an instantaneous basis. To this aim, we have resorted to the following indicators of the rain event: • total rainfall accumulation, For each indicator, we have defined the difference between the daily estimated (denoted with the subscript est) and measured (denoted with the subscript meas) quantity of interest, for each of the 70 considered rain events: Table 5 lists the associated average (on all the events) mean and RMS values of the residuals. The results obtained are satisfactory, as the indicators reveal that the proposed method allows (on average) for estimating the rainfall accumulation with an error in the order of few millimeters. Concerning the detection of peak rain rates, the resulting value of the RMS error is considered to be acceptable, given the difficulty in estimating extreme rain rate values, which might even range between 150 mm/h and 200 mm/h.
Besides the overall statistical indicators, it is worth further investigating the model accuracy by selecting some sample events. Figure 11 depicts the trend of the rain rate during a typical Spring rain event, which is characterized by moderate rain rates and occasional more intense peaks. The equivalent rain height is calculated while using (10) since the event occurred in April (m = 4): the variability of the precipitation is tracked with a significant degree of accuracy, leading to an estimate of the total rainfall accumulation that closely matches the measurements (see Figure 12).         Figures 13 and 14 refer to a typical event in Autumn (November, = 11): in this case, the equivalent rain height is calculated while using (8). The model accuracy is lower than in the previous event, but the results are still quite satisfactory.  The correct evaluation of the equivalent rain height plays a key role in the prediction of the rain rate, as also shown by the curves in Figure 7. Figure 15 depicts the trend of ℎ for the previously mentioned convective and stratiform events, as reported in Figures 11 and 13, respectively, and points out the importance of using local forecast values of ℎ with high temporal resolution. The values in Figure 15 indicate the opposite for the two events, although the rain height is generally higher for convective events than for stratiform ones; in this case, using monthly values of ℎ (or even worst, a single yearly value) would have led to a much lower prediction accuracy than that shown in Figures 11 and 13. The correct evaluation of the equivalent rain height plays a key role in the prediction of the rain rate, as also shown by the curves in Figure 7. Figure 15 depicts the trend of h R for the previously mentioned convective and stratiform events, as reported in Figures 11 and 13, respectively, and points out the importance of using local forecast values of h R with high temporal resolution. The values in Figure 15 indicate the opposite for the two events, although the rain height is generally higher for convective events than for stratiform ones; in this case, using monthly values of h R (or even worst, a single yearly value) would have led to a much lower prediction accuracy than that shown in Figures 11 and 13. points out the importance of using local forecast values of ℎ with high temporal resolution. The values in Figure 15 indicate the opposite for the two events, although the rain height is generally higher for convective events than for stratiform ones; in this case, using monthly values of ℎ (or even worst, a single yearly value) would have led to a much lower prediction accuracy than that shown in Figures 11 and 13.

Comparison with Other Rain Rate Estimation Models: Brief Description
The proposed model is compared with other methods currently available in the literature, which are briefly described hereinafter below.

The NEFOCAST System
NEFOCAST [9] is the name of a project whose purpose is to design a system to detect and monitor rainfall by making opportunistic use of a network of ground terminals: as a matter of fact, such receivers intrinsically have the capability of measuring rain attenuation that affect the satellite downlink, starting from continuous measurements of the received signal-to-noise ratio (SNR). Each

Comparison with Other Rain Rate Estimation Models: Brief Description
The proposed model is compared with other methods currently available in the literature, which are briefly described hereinafter below.

The NEFOCAST System
NEFOCAST [9] is the name of a project whose purpose is to design a system to detect and monitor rainfall by making opportunistic use of a network of ground terminals: as a matter of fact, such receivers intrinsically have the capability of measuring rain attenuation that affect the satellite downlink, starting from continuous measurements of the received signal-to-noise ratio (SNR). Each terminal sends the data to a satellite hub on the return link and such data are then redirected to the NEFOCAST service center, where the rain field maps are generated.
The quantities defined below are expressed as a function of time t. The steps that are required to retrieve the instantaneous rainfall intensity R(t) are the following:

•
obtain the reference SNR level in dry (no-rain) conditions, SNR dry (t); • measure the received signal-to-noise ratio SNR wet (t) and determine whether or not rain is present along the link; • in the presence of rain, extract the corresponding rain attenuation A R (t), that is defined as the decrease in the measured SNR level with respect to SNR dry (t); and, • estimate R(t) from A R (t) as: where a and b are empirical, non-linear regression coefficients (whose values are reported in [9]), θ is the elevation angle of the Earth-satellite link, and h R (t) is the equivalent rain height. The latter quantity is expressed as: where the instantaneous values of h 0 (t) are retrieved from NWP products.

Centralized Method for Real-Time Rainfall Estimation
The model in [8] proposes a real-time and centralized rain rate estimation technique that employs Carrier-to-Noise (C/N) power ratio measurements every five minutes (exploiting both satellite forward-and return-link data), fed back from a network of Ka-band satellite ground terminals to a central gateway, and successively stored in a database. The collected data are then processed while using an Artificial Neural Network (ANN) to detect rain events, and successively rain attenuation is extracted. The authors propose two different solutions to estimate rainfall rate from the inferred rain attenuation: the former takes the Earth-satellite link geometry into account and makes use of the equivalent rain height; the latter is instead a completely empirical relationship. The validity of both approaches is verified while using on-site rain gauges, which report the average rainfall intensity every minute.
The lower case n denotes the n th discrete-time sample of the quantities described below. The estimation procedure can be summarized, as follows: • estimate the baseline C/N level in dry conditions, namely C(n); • detect the presence of a rain event using the ANN described in [8]; and, • once a rain event is detected, determine the corresponding rain attenuation A R (n) as the decrease in the measured C/N (during the rain event) with respect to the last value of C(n) before the beginning of precipitation; • estimate the rainfall rate from A R (n) as: where the values of k and α can be found in [18], θ is the elevation angle of the Earth-satellite link, h s is the altitude of the ground terminal above mean sea level, and h R is the mean annual equivalent rain height, expressed as a function of the mean annual zero-degree isotherm height of the site h 0 : An alternative solution proposed as part of the model is to build an empirical relationship between the observed rain attenuation and rainfall rate, exploiting the training data that were fed to the ANN. The result is an empirical power law of the following kind: where the coefficients a and b are obtained by curve-fitting the used training set to the rain rate measurements produced by the rain gauge.

Model Developed by the Department of Broadband Infocommunications and Electromagnetic
Theory at Budapest University of Technology and Economics (BME-HVT) The model [10] allows for estimating rain intensity from rain attenuation affecting satellite links. No algorithm is provided for retrieving rain attenuation from the received downlink signal: as a matter of fact, the authors have obtained A R by manually identifying the presence of rain events along the satellite link, a process that requires the visual inspection of every received power level time series.
Once rain attenuation is derived, rainfall rate can be estimated by applying the following procedure: 1.
Compute the effective path length L E of the Earth-satellite link, as described in the Recommendation ITU-R P.618 [19]; 2.
Obtain the values of the power-law coefficients k and α, which are available in [18]; 3.
Estimate rain rate as:

Comparison with Other Rain Rate Estimation Models: Accuracy Assessment
This Section presents the comparison of the prediction accuracy that was achieved by the methods described in Section 6.1 and by the proposed model. The following notation will be used to name the methods taken into consideration to simplify the discussion: • NF, which refers to the NEFOCAST system; • CM, defining the centralized method for real-time rainfall estimation [8]. More specifically, the first solution (described by (21)) is labelled as CM1, while the alternative one (described by (23)) is denominated as CM2; and, • BM, associated to the method developed by the BME-HVT. Figure 16 shows the P(R) obtained by using the abovementioned techniques, along with the two curves that were obtained from the measurements and the proposed model, respectively. In addition, the corresponding statistical indicators are reported in Table 6. 1. Compute the effective path length of the Earth-satellite link, as described in the Recommendation ITU-R P.618 [19]; 2. Obtain the values of the power-law coefficients and , which are available in [18]; 3. Estimate rain rate as: (24)

Comparison with Other Rain Rate Estimation Models: Accuracy Assessment
This Section presents the comparison of the prediction accuracy that was achieved by the methods described in Section 6.1 and by the proposed model. The following notation will be used to name the methods taken into consideration to simplify the discussion: • , which refers to the NEFOCAST system; • , defining the centralized method for real-time rainfall estimation [8]. More specifically, the first solution (described by (21)) is labelled as 1, while the alternative one (described by (23)) is denominated as 2; and, • , associated to the method developed by the BME-HVT. Figure 16 shows the P(R) obtained by using the abovementioned techniques, along with the two curves that were obtained from the measurements and the proposed model, respectively. In addition, the corresponding statistical indicators are reported in Table 6.   The results indicate that the proposed model shows the highest accuracy. A critical aspect that is common to all of the models listed above is the limitation in estimating high rain rates: all the curves diverge significantly from the reference one (red line), as the rain rate increases. This issue is significantly lower for the proposed model, which introduces a more complex and accurate method for estimating the equivalent rain height. The lower accuracy of the other models is likely due to the following reasons: • NF does not consider any path reduction factors, which, on the contrary, are included in the prediction model in [19], and are intended to model the spatial inhomogeneity of rainfall along the path.
• CM1 presents similar problems as NF; in addition, the equation γ R = kR α is inverted while using the power law coefficients provided by the recommendation ITU-R P.838-3, which differ from the one locally derived for Milan (see Figure 6); • The limited accuracy of CM2 is mainly ascribed to its strongly empirical nature; • BM reveals the same issues as CM1. In addition, the inversion procedure requires as input information on the long-term meteorological statistics of the site (such as R 0.01 ).
Tables 7-9 complete the model performance evaluation by reporting the statistical indicators for R T , R M , and R. The results confirm the higher accuracy of the proposed model.

Conclusions
This work aimed at developing a rainfall rate estimation algorithm, being composed of two main parts: a rain attenuation to rain rate conversion function, which receives the instantaneous rain attenuation and equivalent rain height as input, and provides the estimated time series of the rain rate as output; a methodology for estimating the equivalent rain height h R , based on the 0 • isotherm height and on the local monthly trend of the convectivity index.
The model was applied to the Ka-band propagation data that were collected in Milan in the framework of the Aldo Paraboni experiment, and its performance was assessed by using the rain rate data collected by the collocated disdrometer as reference. The performance evaluation was carried out both on a statistical basis and on an instantaneous basis, the latter consisting in the inspection of the individual time series. If compared to other methodologies, the proposed model achieves the best performance in estimating the rain rate statistics, with an overall root mean square of the prediction error equal to 1.19 mm/h. The additional error indicators calculated for individual events confirms the model performance: in all cases the proposed model shows to the best prediction accuracy, with RMS values equal to 5.34 mm, 11.83 mm/h, and 1.52 mm/h when the rainfall accumulation, the peak rain rate, and the average rain rate are considered, respectively.
Notwithstanding the very good results achieved, it is worth mentioning that the proposed model was calibrated and tested on the same reference dataset: additional data, not yet currently available to the authors, are required for corroborating the model applicability and performance in other sites.