Retrieval of rainfall fields in urban areas using attenuation measurements from mobile phone networks : A modeling feasibility study

Rainfall monitoring is an important global issue in urban hydrological applications, such as flood warning and water resources management systems. Until the present time, rain gauges and weather radars have been widely used as sensors to provide rainfall information with a detailed resolution; most cities in the world however are inadequately equipped. Recently, commercial microwave links (MWL) have been proposed as a new means of monitoring space-time rainfall. A 15 transmitted signal along such links is known to be attenuated by rainfall, hence the measurement of this signal attenuation could serve to estimate path-averaged rainfall intensity. The density of commercial MWL is typically high in most cities today, which raises new questions over the possibility of retrieving rainfall using signal attenuation data from multiple links. The objective of this article is to assess the feasibility of retrieving rainfall fields in urban areas using rain attenuation data from commercial MWL that are mainly operated by mobile phone companies. This work is based on a simulation framework 20 applied to a real case study. The study area is the city of Nantes, France. Rainfall datasets containing 207 weather radar images recorded by the Météo-France Agency's C-band at high spatial (250 m x 250 m) and temporal (5 min) resolutions are first used to generate rain attenuation data over the existing mobile phone network, which combines 256 microwave links operating at 18, 23 and 38 GHz. The rain attenuation data generated are used as a real signal dataset. A novel retrieval algorithm is then proposed to convert the rain-induced attenuation data into a rainfall map. A priori knowledge introduced to 25 initialize the algorithm heavily influences retrieval performance if the problem to be solved is under-determined, as is the case herein. The capabilities as well as limitations of the retrieval algorithm, as regards capturing different rainfall variability, are evaluated. A detailed sensitivity analysis, carried out with respect to various parameters including a priori knowledge, decorrelation distance, and the retrieval performance of the algorithm depending on the density level of the MWL network is 30 also evaluated in a light rain, a shower and amidst storm events. The conclusion, based on 200+ retrieval tests, states that the Hydrol. Earth Syst. Sci. Discuss., doi:10.5194/hess-2016-540, 2016 Manuscript under review for journal Hydrol. Earth Syst. Sci. Published: 21 October 2016 c © Author(s) 2016. CC-BY 3.0 License.


Introduction
Rainfall may constitute an obstacle to radio wave propagation in the atmosphere at frequencies ranging from 5 to 50 GHz (Hogg, 1968;Atlas & Ulbrich, 1977).This inconvenience offers new opportunities to measure rainfall.The measurement of rainfall rates using radio links set up by phone companies has become an active research topic.These links mainly operate at frequencies attenuated by rainfall; moreover, measuring the path-integrated attenuation between the link transmitter and receiver could yield the mean rainfall rate along this path, as suggested by (Upton, et al., 2005;Messer, et al., 2006;Leijnse, et al., 2007).Several experiments, conducted in very different contexts, have confirmed these results: (Doumounia, et al., 2014) obtained encouraging results in Sahelian West Africa along a 29-km link at 7 GHz; (Fenicia, et al., 2012) tested different analysis methods on two dual-frequency links in Luxemburg City; (Ostrometzky, et al., 2015) proposed a method that allows estimating all types of precipitation quantities.Recent experiments, involving larger datasets recorded on a greater number of telecommunication links and validation by rainfall measurement devices (rain gauges, disdrometers), have provided robust assessments of these initial findings (Overeem, et al., 2011;Rayitsfeld, et al., 2012).These promising results must not neglect however the presence of several sources of errors, whose influence has been studied by the following authors: small-scale variability of rainfall along the links (Berne & Uijlenhoet, 2007), attenuation due to wet antennae (Leijnse, et al., 2007), spatial variations in DSD (Drop Size Distribution) along the link (Leijnse, et al., 2010;Zinevich, et al., 2010), all of which present a very comprehensive study of the various sources of error and uncertainty.
Telecommunication links might also be used to improve the operations of conventional rainfall measurement devices: detection of malfunctioning rain gauges in urban areas (Bianchi, et al., 2013), determination of the mean field adjustment of weather radar data (Cummings, et al., 2009), and the detection and correction of attenuated radar data (Rahimi, et al., 2006).
One important and interesting aspect is the mapping of rainfall fields from path-integrated attenuation measured on a network of telecommunication links.This subject was first addressed by (Giuli, et al., 1991), who proposed a "microwave tomographic inversion technique" that considers rainfall fields as a linear combination of basis functions.The authors also tested the method within the academic framework of simulated rainfall fields sampled by a regular network of radio links.
Their method has since been further improved: (Giuli, et al., 1999;Cuccoli, et al., 2011) proposed reconstructing rainfall fields from path-integrated attenuation measurements on radio communication links by means of a tomographic processing Hydrol.Earth Syst.Sci.Discuss., doi:10.5194/hess-2016-540,2016 Manuscript under review for journal Hydrol.Earth Syst.Sci.Published: 21 October 2016 c Author(s) 2016.CC-BY 3.0 License.method referred to as the "Combined Deterministic-Stochastic Retrieval Technique".These authors obtained valuable results when testing their method using a simulation approach applied to a real case study, namely rainfall fields observed by radar and the actual linked network of the city of Florence (Italy).(Zinevich, et al., 2008) used the SIRT, i.e. the simultaneous iterative reconstruction technique described by (Kak & Slaney, 1988), to reconstruct rainfall fields in a study area equipped without spatial regularity or frequency homogeneity.(Zinevich, et al., 2009) combined a space-time advection model of rainfall with an extended Kalman filter (Welch & Bishop, 2006) to perform the retrieval from path-integrated attenuation measurements.(Goldshtein, et al., 2009) proposed a stochastic interpolation algorithm based on the Inverse Distance Weighting function (Shepard, 1968) to reconstruct rainfall fields over a sparse and arbitrary cellular network geometry.Most of the cited studies were intended to retrieve rainfall in urbanized areas well equipped with telecommunications links.On a very different scale, (Overeem, et al., 2013) implemented an extensive network of 2,400 links to retrieve the space-time dynamics of rainfall for the 35,500  2 of the Netherlands, by interpolating the rainfall estimates derived at link locations using the Kriging method (Creutin & Obled, 1982).(Gosset, et al., 2015) emphasized the potential benefit of commercial MWL for mapping rainfall in areas with poor or nonexistent coverage offered by conventional measurement devices (e.g.rain gauge, weather radar), including in urban areas.This article addresses the special case of cities well equipped with MWL operated by telecommunications companies.
A small number of studies have been devoted to reconstructing rainfall fields from path-integrated attenuation on these links; the feasibility of this reconstruction process has not yet been fully established.In the follow-up to (Cuccoli, et al., 2011), the objective has been to assess the feasibility of estimating rainfall fields from attenuation measurements on microwave communication links by means of a simulation approach based on a real-world case study.The series of observed weather radar images, recorded at a very high spatial (250 m) and temporal (several minutes) resolution, are used to simulate attenuation measurements along the communication links of the Nantes metropolitan area.These attenuation "data" will then be applied to retrieve the rainfall field by running an inverse algorithm.The retrieval efficiency will be evaluated by comparing the observed and retrieved rain fields.A highly detailed analysis of sensitivity to influential factors (i.e.available data and associated errors, inverse algorithm, a priori information, spatial and temporal rain field structure, communication link network geometry) will be performed in order to assess the feasibility of the proposed method.Section 2 will present the case study under consideration as well as the simulation of path-integrated attenuation measurements.Section 3 will then discuss the method used to retrieve rainfall fields and its application conditions in the case of the studied problem.The results obtained and a sensitivity analysis conducted on these results will be addressed in Section 4.

Methodology description
This section will present both the problem raised by the retrieval of rain rate maps from path-attenuated measurements on an MWL network covering an urban area extensively.The corresponding attenuation model will be formulated in Subsection 2.1, while Subsection 2.2 will examine the inverse algorithm used to solve this retrieval problem, by relying on the works of (Menke, 1989) and (Tarantola & Valette, 1982a).Subsection 2.3 will present the conditions for applying this algorithm to the problem at hand.

Path-integrated attenuation by a network of microwave links
Microwave attenuation by rainfall along a link path can be expressed as: with A l being the path-integrated attenuation (hereafter denoted PIA) along link l, (T l , R l ) the transmitter and receiver coordinates of this link, R(u) the rainfall intensity at location u, and (a l , b l ) the parameters of an attenuation law dependent on the link wavelength, polarization, temperature and ultimately on the DSD of rainfall (Olsen, et al., 1978).
Let the rainfall field be defined as a pixelated image, e.g. a weather radar image, then Eq. ( 1) can be expressed as follows: where n(l) is the number of pixels in the rain image crossed by link l, L k(l,j) the path length of link l in pixel k(l, j) of the rain image, and R k(j) the rainfall intensity in pixel k(j, l).The function k(l, j) serves to associate the rain image pixels with the MWL that cross them.
Considering that the MWL network contains N links, the overall model relating the PIA on all the links to the pixelated rain image can be written as: where the vector combining the rainfall rate on the M rain field pixels crossed by a link, and ε = [ε 1 … … … ε N ] the vector combining the modeling errors on all N links.Let's note that the rain pixels involved in model m do not comprise all rain image pixels.The number of uncrossed pixels depends on the microwave link characteristics and the spatial resolution, ∆ = ∆x = ∆y , of the rain pixels, which are assumed to be square.

The identification algorithm
The objective of this algorithm is to retrieve the rainfall field that best reconstitutes the measured PIA values of the MWL network according to the model defined by Eq. ( 4).The method chosen to solve this inverse problem is the conventional inverse approach, as detailed in (Menke, 1989), while the algorithm used in this work is the one proposed by (Tarantola & Valette, 1982a).The solution minimizes the following expression: where Φ is the likelihood function, t signifies transpose,  0 denotes the vector of measured PIA values,   is the a priori rainfall rate on the pixels crossed by a microwave link, and   and   are the covariance matrices of residuals between truth and a priori values of  and , respectively.The statistical distributions of both [ −  0 ] and [ −   ] are assumed to be unbiased and Gaussian.(Menke, 1989) demonstrated that the solution vector ′ satisfies: Where  is the matrix of (first-order) partial derivatives of model m.If model m is nonlinear, Eq. ( 5) can then be solved using an iterative method, i.e.: in which   constitutes the result of the k th iteration and   the matrix of partial derivatives of model m at point   .Further information about the stability, convergence and uniqueness of the solution of such nonlinear problems can be found in Chapter 9 of (Menke, 1989) and in Tarantola and Valette (1982a and b).The general case of inverse problems involving non-Gaussian statistics has been addressed by (Tarentola, 1987).

Remarks
Before applying this algorithm to the retrieval of rainfall fields using PIA measurements from an MWL network, let's consider the a priori information and nature of the studied problem.The a priori information is composed of: i) the vector   and its covariance matrix   , which characterize the level of confidence in the data; and ii) the vector   and its covariance matrix   , which contain initial knowledge of the parameters to be identified (i.e.rainfall field) and the level of confidence ascribed to this knowledge.The solution output by the inverse algorithm results from a compromise between two extreme states: i) a solution that perfectly fits the observed data through the theoretical model; and ii) a solution remaining very close to the a priori information on the parameters.The assistance provided by this a priori information actually depends on the nature of the problem.If the observed data are insufficient or if confidence in the data is weak, the problem is then considered underdetermined and the a priori information assumes a dominant role.On the other hand, if the problem is over-determined (i.e.very high quality data in sufficient quantity), the solution no longer depends on the a priori information.The special case of rainfall field retrieval from PIA measurements in an MWL network depends on several factors: -The MWL network heterogeneously covers the study area, as illustrated in Figure 1, which shows that the link density decreases from the city center out to the suburbs; -The spatial resolution adopted to represent the rainfall field complicates to some extent problem resolution.
Depending on this spatial information, the number of pixels uncrossed by one or more links may vary substantially.
Yet the increased pixel size can lead to a loss of representativeness if the rainfall field displays strong spatial variability.In summary, it appears that the problem to solve is "heterogeneously" determined and that both practical implementation of the algorithm and choice of the a priori information are capable of significantly influencing the solution.

Application of the retrieval algorithm
This section will present the conditions for applying the algorithm to retrieve the rainfall field in an urban area through use of a dense MWL network.

Attenuation data and associated covariance matrix
The errors associated with the attenuation measurement along HF links and the rainfall rate estimation from these measurements have been modeled by (Zinevich, et al., 2010) and studied in detail by (Leijnse, et al., 2010).The vector   encompasses PIA measurements on the N links.The attenuation measurements are assumed to be mutually independent.
Consequently, the covariance matrix   is assumed to be diagonal, with its diagonal term being expressed as: with the  parameter being equivalent to the coefficient of variation of the statistical distribution of the assumed measurement error (Zinevich, et al., 2010) on the PIA values; the 0.01 constant has been added to prevent the matrix from becoming singular.

A priori values of rainfall
The a priori knowledge on model parameters is described by the rainfall rate vector   and its covariance matrix   .The selected approach consists of deducing the a priori rainfall at rain pixel i from the measured PIA values on the links closest to this location, i.e.: where   is the number of considered links;  0 () the measured PIA on link (), which is one of the   links closest to pixel ; and  () and  () are parameters of an R-A relationship adapted to link ().The distance between a link and a pixel is derived from the orthogonal projection of the pixel (i.e.center) on the link.
The a priori error variance of  0  is estimated from the variance of the   values of   used to estimate  0  ; it is expressed as: and the terms of the a priori covariance matrix   are estimated by: with   being the distance between rain pixels i and j, and the D parameter controlling the decorrelation between errors on the a priori rain rates.
It is observed that the number   of the closest selected links assumes significant importance during algorithm initialization.
Let's consider that at   = 5, the initialization remains very local.The maximum value of   equals the total number of study area links, hence each rain pixel is assigned the same a priori value (and error variance), which is equal to the mean value and the variance of the rainfall field derived from the measured PIAs.Moreover, let's note that the initialization is not homogeneous over the entire study area since the neighborhood corresponding to the   closest links depends on the link density.

Application of the algorithm -Grid nesting procedure
Subsection 2.3 pointed out that available information decreases from the city center to the suburbs.In order to take into account this specific feature, it is proposed herein to apply the identification algorithm under the following conditions: i) the pixels of the rain field targeted by the identification are merely those crossed by MWL; and ii) the pixel area is adapted to the MWL density.
The identification algorithm is then implemented according to the following grid nesting procedure: Step 1: Rain field retrieval is performed first at a low spatial resolution (22  2 ) over the largest possible geographic area covered by the link network; Step 2: The solution obtained serves as a priori for the rain field retrieval at an improved spatial resolution (11  2 ) over a subarea limited to the pixels crossed by the links at this level of resolution.During this step, the uncrossed pixels retain the values assigned in the previous step and remain a solution; Step 3: This procedure can ultimately be repeated to further improve the spatial resolution in once again restricting the subarea to the remaining crossed pixels.At this point in our case, the spatial resolution is set at 0.50.5  2 , which is acceptable as a final spatial resolution of the retrieval process.
Practically speaking, this iterative grid nesting procedure calls for refining rainfall field retrieval as the link density increases.The iterative procedure explained above is illustrated in Figure 2.

Presentation of the case study: City of Nantes
This study seeks to assess the feasibility of retrieving rainfall fields in urban areas from existing MWL deployed by telecommunications companies.Our case study is based on the situation in the city of Nantes.

Microwave links
The city of Nantes, France is covered with a set of MWL owned by the four mobile communications companies operating in France.At the end of 2012, the urbanized part of this city was equipped with a total of 256 MWL (www.cartoradio.fr)operating at 18, 23 and 38 GHz (as shown in Fig. 1).The lengths of these microwave links range from 0.3 to 16.8 km, with an average of 4.2 km.The main characteristics of this set of links are detailed in Table 1.The rows labeled "count", "mean", "min" and "max" list respectively the number of microwave links, their average, minimum and maximum lengths.Within the network service area, nearly 95% of the MWL at 38 GHz are shorter than 5 km.This same percentage threshold reveals almost 11 km and 10 km for the links at 18 and 23 GHz, respectively.These data provide a clear indication of the fact that microwave link lengths at higher frequencies are basically assumed to be shorter and vice versa.At all frequencies, the minimum lengths are practically the same, i.e. at less than 1 km.

Classification of microwave link density
It is important to define the density level of the microwave link network, as this parameter allows differentiating between dense and sparse regions within the study area.The reason for such a distinction is to understand the influence of network topology on rainfall retrieval accuracy.The density level of the MWL network has therefore been computed.The result of this computation will be called a "pixel density map", which can then be used to validate network system performance across various regions involved in the rainfall retrieval process.
A basic assumption herein is that the higher link density should yield more weight on the surface.The following expression is introduced to obtain the  th pixel area density: where  is the link index,  the pixel index,  the number of links, W the pixel density vector,  the length of the intersected part of link  (in km).If Eq. ( 11) were applied to all pixels crossed by microwave links, the pixel density vector would be obtained.This vector could then be reshaped into a two-dimensional pixel density map.
Figure 3 shows a histogram of the pixel density computed over 256 MWL within the study area.The distribution of W in this figure is heavily right-skewed, indicating that the average value of W is greater than its median.The density level is considered "low" if the condition  <   is satisfied; in our case,   = 0.5.Similarly, the density level is considered "moderate" if W lies between 0.5 and 1.Beyond the moderate level, density is considered sufficiently "high" since it is assumed to be acceptable if some 30% of W values exceed 1.However, this percentage depends on the user's choice.Figure 4 shows the pixel density map derived at a 22  2 resolution (i.e.pixel size) over the entire MWL network area.Note that the density levels are classified in the colors gray, red and blue, which represent low, moderate and high density regions of the MWL network, respectively.sufficient to cover the area of microwave links.For purposes of illustration, the radar rain maps for each considered event are displayed in Figure 5.
A total of 34 rain periods recorded between 2009 and 2012 have been selected for this study.According to (Emmanuel, et al., 2012), the days of recordings were partitioned into four groups of different rainfall field types.Thanks to a meteorological analysis combined with radar image visualization, three categories of meteorological situations could be detected at the mesoscale: warm sectors, fronts, and the tail end of low pressure systems.These situations yield four distinct types of rainfall fields, namely: 1) light rain (correlated with warm sectors), 2) storms organized in rain bands (correlated with fronts), 3) storms less well organized (correlated with fronts), and 4) showers (correlated with the tail end of low pressure systems).(Emmanuel, et al., 2012) studied the spatial variability of these 4 types of rainfall, which display distinct characteristics in terms of intermittency and rain variability.In our study, 7 rain events were extracted from these periods by collating a total of 207 rainfall maps, representing these 4 rainfall types characterized by different spatial and temporal variabilities.Descriptions of the duration, period and number of rainfall maps for each event are provided in Table 2.

Path-Integrated Attenuation data
To generate PIA data over the given MWL network, the following steps have been performed: -Extract a rainfall map of weather radar and superimpose it on the study area.In this case, the pixel index of the rainfall map is expected to be the same as the index of the discretized area; -Compute the PIA along each link using Eq.The output of the above steps is to be substituted for the real rain attenuation data potentially obtained from telecommunications companies.The sources of error during the rain attenuation measurement process are twofold: i) the PIA model (Eq.3); and ii) instrumental impairment such as quantization, antenna wetting and temperature, and gas effects.Since quantifying the error source with precision lies outside the scope of this paper, a simple approach is considered for simulating the nature of the real signal attenuation caused by rain.This approach has been introduced using a Gaussian distributed random variable with variance equal to 5% of PIA attenuation, as computed along each link of the MWL network.(Atlas & Ulbrich, 1977) showed that uncertainties due to drop size distribution variability when determining the  (Atlas & Ulbrich, 1977).Higher percentages indicate that our confidence in the accuracy of rain attenuation measurement is lower.However, a higher percentage (i.e.20%) has also been tested to show that the attenuation measurement might be influenced by greater sources of error as well.Moreover, the signal quantization effects are introduced after computing rain attenuation in the presence of the PIA model-related error source.This quantization error occurs when the signal is truncated at a certain level of precision.Typical precision values are 0.1 and 1 dB, at which the actual signal becomes truncated.The signal quantization error variance remains uniform and equal to σ q 2 = ∆ 2 12 (Widrow & Kollàr, 2008).The magnitude of the quantization effect has been well understood in a number of studies (Goldshtein, et al., 2009;Zinevich, et al., 2010), with this type of error source also exerting a significant influence on shorter links during lower-intensity rain events.The higher quantization value indicates less accuracy and precision inherent in the antenna station equipment.
The steps described above are applied to generate the rain attenuation data and then repeated for all 207 rainfall maps.Figure 6 illustrates the attenuation dataset produced in light rain, showers, organized and unorganized storm events.The y axis represents the average path-integrated attenuation values along each link derived on the 113, 26, 34 and 34 rainfall maps associated to those events, respectively.Here again, the data generated have been classified according to link length in 3 different ranges, namely (0, 2], (2, 5] and (5,17] km, which practically represent the percentile intervals of (0, 25%], (25%, 50%] and (50%, 75%] of all MWL lengths, respectively.This classification is useful in providing a deeper understanding of the rain attenuation variability derived depends on link length from one rain event to the next.
It is evident that the PIA value gradually increases with increasing link length and frequency.Another key finding is the rising trajectory of this value at higher rainfall variability, especially in the unorganized storm case (i.e.lower right figure), which reveals a "close-to-reality" type of measurement error sources.In this study, the simulated error sources added to the PIA measurement process generally reflect the nature of possible error sources discussed in experimental and theoretical studies (

Sensitivity analysis
The sensitivity analysis offers insight into understanding the interrelationships between algorithm parameters.The significance of algorithm sensitivity is to test the influence of various parameter values on the solution vector.In other words, the result of this sensitivity test reveals that the algorithm is either sensitive or insensitive to changes in algorithm parameter values.
In this section, the analysis of sensitivity to algorithm parameter values will be discussed.Subsection 5.1 will describe the sensitivity analysis protocol, while Subsection 5.2 will address the influence of the a priori parameters (i.e.  and D) on rainfall retrieval.

Sensitivity analysis protocol
The overall structure of the sensitivity analysis testing the algorithm parameters is given in Figure 7. Levels 1 and 2 denote the values of sensitivity parameters to be tested, which consist of a priori rainfall selection   and decorrelation distance D, respectively.The Level 3 procedure is performed in the following order: -Rainfall map: Extract the rainfall map from the weather radar.Map resolution is 250 m x 250 m, as obtained from the radar data (see Section 4); -Rain attenuation data generation: Rain attenuation data over the microwave link network are generated according to the procedure set forth in Subsection 4.3; -Rainfall retrieval: The grid nesting procedure described in Subsection 3.1 is applied to perform rainfall retrieval; -Comparison: The rainfall map retrieved by the proposed algorithm is compared with an observed radar rainfall map by employing the NSE criterion commonly used in the field of hydrology (Nash & Sutcliffe, 1970), expressed as follows: The NSE criterion provides an overall assessment by taking into account both the observed (  ) and retrieved (  ) rainfall maps, with n being the number of pixels in the compared maps.The interval for NSE values lies between −∞ and 1: the closer the NSE value is to 1, the more accurate the model retrieval.
Level 3 is applied to the 207 rainfall maps, which represent light rain, showers, organized and unorganized storm events.
Let's note that each event has been tested and analyzed separately.This protocol helps determine the sensitivity of retrieval algorithm parameters to various types of rainfall.

Sensitivity to the a priori rainfall map
The objective here is to test the influence of the a priori rainfall.Each rain pixel is set at a value equal to the mean of the rainfall estimates, as derived from PIA values on the links closest to its location, i.e.   .These sensitivity tests pertain to the neighborhood with   spanning a range that includes 5, 7, 10, 15, 100 and 256, i.e. from a very local to the global neighborhood, respectively.Let's recall that the standard deviation in a priori variance associated with the a priori rainfall has been set at the variance of the considered rainfall estimates.Applying the sensitivity protocol described in Subsection 5.1 yields Figure 8, which illustrates the results of sensitivity to the a priori parameters.The red dots denote the highest NSE values, which indicate that the algorithm performs best with these choices.Increasing the values of   certainly decreases algorithm accuracy regardless of rainfall variability.One important feature found in these results however is that higher rainfall variability widens the optimum   parameter value ranges (see, for example, the difference in the number of red dots between light rain and unorganized storm events).

Sensitivity to the decorrelation distance of the a priori rainfall error
The purpose here is to test the sensitivity of the retrieval algorithm to parameter D, see Eq. (10).The results of the analysis in Figure 8 show that the algorithm is also sensitive to the value of D: it appears that the influence of D increases with rainfall variability.As opposed to the   parameter, the optimal D value is found to lie in a tight range, i.e. very few choices, regardless of rainfall type.

Evaluation
This section will evaluate the capabilities and limitations of the proposed rainfall retrieval algorithm.A series of rainfall events recorded by weather radar are used as a basis to assess the proposed algorithm.This evaluation is performed at a 0.5 km x 0.5 km resolution with respect to the optimum (selected) parameters of the algorithm given in Table 3.The evaluation principle and conditions will first be described in Subsection 6.1; Subsection 6.2 will then discuss the results obtained.

Evaluation principle
The evaluation principle is defined as follows: i) an efficiency curve is calculated for each event (rainfall type) that allows assessing the retrieval performance of the algorithm vs. rainfall variability; and ii) the evaluation zone is classified in order to verify the retrieval efficiency in different regions of the MWL network (this process is performed based on the pixel density map given in Subsection 4.2).
The initial findings of this evaluation protocol consist of: a) an overall comparison between the retrieved (by the algorithm) and the observed (by weather radar) rainfall maps; and b) a separate analysis of the influence of the inverse algorithm and a priori knowledge on retrieval performance.These features help understand the added value of a priori knowledge depending on rainfall variability.Subsequent findings stem from evaluating the method over the entire monitoring network system and evaluating retrieval efficiency on the basis of microwave link network density.This latter evaluation relies on the pixel density map (as detailed in Subsection 4.2), which depicts the microwave link density level.

Evaluation results
Regarding the initial findings of our analysis, examples of retrieved rainfall for light rain, showers, organized and unorganized storm events are shown in Figure 9.These maps were obtained using Level 3, which is the retrieval protocol depicted in Figure 7.It is evident that the spatial variability for all rainfall types has been far better captured, especially for the case depicted in the second row of this figure.Since the optimized algorithm parameters resulting from the sensitivity analysis have been used in these tests, the a priori initialization procedure heavily influences the algorithm output.It is apparent however that the algorithm smoothes the solution in all cases, as compared to the observed maps.
Local initialization has been shown to significantly improve retrieval algorithm performance.Yet, two critical questions arise: i) does this improvement stem solely from the a priori knowledge; and ii) does applying the inverse algorithm benefit rainfall retrieval.To answer these questions, we have compared the a priori knowledge with the inverse algorithm performance at each step of the grid nesting procedure (Fig. 2) individually, meaning that the rainfall map values retrieved by the algorithm have been compared with the a priori rainfall map values introduced before the actual retrieval.This comparison was performed at each step of grid nesting.Figure 10 illustrates the results of the comparison using the NSE criterion, as plotted in efficiency curve.Let's note that the figures in the 1 st , 2 nd and 3 rd columns result from the 1 st , 2 nd and 3 rd steps of the grid nesting procedure, respectively.It is observed that the added value of the algorithm is highly significant at 22  2 , especially in the presence of strong rainfall variability, see for example the shower event.In other words, the greater the rainfall variability, the larger the influence of the inverse algorithm on output.This influence then moderates and becomes less significant at resolutions of 11  2 and 0.50.5  2 , respectively.Interestingly, the overall contribution of the inverse algorithm becomes roughly the same as the retrieval resolution decreases from 22  2 to 0.50.5  2 .This outcome is due to the fact that the a priori solutions in the 2 nd and 3 rd columns do not actually represent the local initialization.Instead, they are the results of the algorithm run during previous steps (i.e.steps 1 and 2 in the grid nesting procedure) and accepted as the a priori rainfall for the subsequent step.
Figure 11 displays the overall retrieval efficiency for the algorithm in NSE criterion vs. microwave link density.This plot is derived from the pixel density map.It is important to note that the line colors in Figure 11 are the same as those shown in Figure 4, except for the yellow color, which has been labeled as "whole density" (thus indicating an evaluation result without the density map).It would seem that the algorithm accuracy within the denser part of the network, i.e. blue line, is quite high compared to the sparser part, i.e. the rainfall retrieved in the city center is more accurate than that in the suburban zone.
Regions with a moderate density label (red line) however indicate an unusual performance in the first few retrieved maps.The explanation for such behavior is that the number of retrieved pixels to be compared with observed ones is much less than the pixels elsewhere, i.e. the red or blue regions.
The scatter plots of the retrieved rainfall shown in Figure 12 illustrate another perspective on the performance of the proposed algorithm.Regression lines in black and red colors have been placed on the observed maps (by weather radar) and on the estimated maps (by the algorithm), respectively.A systematic underestimation is observed across all events for both the whole and high density regions within the study area.However, this underestimation is decreasing in the denser regions of the network, as seen in the figures in the right-hand column.Note that higher rainfall rates are always being smoothed, thus resulting in an underestimation.This smoothing process is explained in two ways: i) the solution vector retrieved during the 1 st step of the grid nesting procedure essentially belongs to rural areas (moreover, it is not subsequently updated by the algorithm in the next steps); and ii) the longer microwave links are also causing an underestimation.

Conclusion
The objective of this study has been to assess the feasibility of rainfall mapping by means of an inverse method applied to signal attenuation data generated on an existing microwave link network within a simulation framework.The new retrieval algorithm based on the inverse method has been proposed to demonstrate the capability of microwave links for rainfall monitoring in urban areas.We have demonstrated that the a priori knowledge used to initialize the algorithm heavily influences retrieval accuracy even when rainfall variability is high.A detailed sensitivity analysis focusing on the a priori parameters has also been performed to test algorithm stability at different levels of rainfall variability.The algorithm applied has been found stable enough even if the problem to be solved remains overly under-determined, as found to be the case herein.Various evaluation procedures have been carried out to analyze monitoring system performance.One weakness of the algorithm identified during the evaluation process is its systematic underestimation, especially with very high rainfall intensity.However, the proposed algorithm is still suitable with larger-magnitude measurement errors.

Future work
This research study opens several new avenues of future study to better understand the rainfall retrieval capability of commercial microwave links.
Further work could be dedicated to assessing the long-term effects of the cellular network used for rainfall monitoring based on a more extensive dataset, e.g. one week of instantaneous rainfall provided by weather radar.In addition, the structure used to model uncertainty introduced into the algorithm could be changed based on the same recommendation given in the previous section, in order to provide a deeper understanding of algorithm capabilities as well as limitations.Since this retrieval algorithm relies heavily on a priori knowledge, different types of a priori initialization procedures can as an option be introduced to raise the degree of accuracy in the output.As regards the simulation framework adopted for this work, further research is required to improve the quality of the resulting rain attenuation data in order to reflect a real type of signal data capable of being measured at microwave antenna stations.Moreover, the error magnitude, i.e. standard deviation, introduced in the algorithm does not take modeling errors into account.This shortcoming may be remedied however by devising a new structure that includes link characteristics (essentially frequency and length) and climate conditions, such as temperature, humidity, gas and water vapor effects.Such an error structure could be developed on the basis of the methodology for modeling error sources and uncertainty proposed by (Zinevich, et al., 2010).
Furthermore, it would be interesting to test these two proposed retrieval algorithms in different cities throughout the world based on an actual case study that requires signal attenuation data obtained directly from cellular phone companies.This idea however will naturally raise the challenge of having to cleanse the signal data recorded at antenna stations before use in rainfall retrieval.Signal attenuation unrelated to rainfall includes baseline effects (i.e.gas, water vapor, temperature, wind) and antenna wetting following the rainy period.The algorithms developed by ( A global perspective on the main objective of this research might entail combining the commercial microwave links with weather radar in order to improve the accuracy of real-time rainfall monitoring within urban areas. Hydrol.Earth Syst.Sci.Discuss., doi:10.5194/hess-2016-540,2016 Manuscript under review for journal Hydrol.Earth Syst.Sci.Published: 21 October 2016 c Author(s) 2016.CC-BY 3.0 License.

4 . 3
Hydrol.Earth Syst.Sci.Discuss., doi:10.5194/hess-2016-540,2016 Manuscript under review for journal Hydrol.Earth Syst.Sci.Published: 21 October 2016 c Author(s) 2016.CC-BY 3.0 License.Rainfall data Rainfall data have been provided by the C-band weather radar at Treillières, with a 3-dB beam width equal to 1.25 deg.This radar station, located 10 km north of Nantes, is labeled as "C-Band" in Figure 1.Radar images are available at 5-min intervals, with an area size of a 128 km x 128 km Cartesian grid at a spatial resolution of 250 m x 250 m.Since the radar image size (i.e.128 km x 128 km) is larger than the study region, only a 40 km x 40 km piece has been selected, which is

( 3 )
in the presence of measurement error ε.The a and b coefficients of the PIA model in Eq. (3) at 18, 23 and 38 GHz are obtained from (ITU-R, 2005); -Quantize the computed PIA at a 0.1-dB resolution.
Hydrol.Earth Syst.Sci.Discuss., doi:10.5194/hess-2016-540,2016 Manuscript under review for journal Hydrol.Earth Syst.Sci.Published: 21 October 2016 c Author(s) 2016.CC-BY 3.0 License.path-averaged rain rate increased from 10.8% to 39.7% as frequency decreased from 34.86 to 9.37 GHz.These authors found however that should the k-R relation be quasi-linear, i.e. b ≈ 1, then this influence would be negligible, resulting in less than 10% average error.According to the authors, the k-R relation linearity reaches a peak at 35 GHz.This finding is clearly applicable to our case since the exponent values (i.e.b coefficient) in Eq. (3) at 18, 23 and 38 GHz are nearly 1.More recently, theoretical as well as experimental studies (Leijnse, et al., 2007; Leijnse, et al., 2010; Zinevich, et al., 2010) have also confirmed the findings reported by Hydrol.Earth Syst.Sci.Discuss., doi:10.5194/hess-2016-540,2016 Manuscript under review for journal Hydrol.Earth Syst.Sci.Published: 21 October 2016 c Author(s) 2016.CC-BY 3.0 License.

Figure 1 .
Figure 1.Mobile phone network located in Nantes city.The zoomed-in part of the figure (yellow square) represents the location of microwave link antennas (Transmitter and Receiver stations) operating at 18, 23 and 38 GHz.Two yellow arrows are also depicted as an example to indicate transmitter stations.10

Figure 2 :
Figure 2: Grid nesting procedure applied to one pixel crossed by a link.The left and the right columns represent the a priori and retrieved rainfall by the algorithm, respectively.The rows (top-down direction) are the stage 1, 2 and 3 in grid nesting procedure.R stands for rain rate value.The segment line (in blue) is a sample microwave link and its labels, `Tr.' and `Rec.',stand for 10

Figure 3 .
Figure 3. Histogram of the pixel density map values.

Figure 4 .
Figure 4. Pixel density map computed over the microwave link area.

Figure 6 .
Figure 6.Boxplot of overall generated PIA data set.Note that y axis represents the mean values of path integrated attenuation.The figures (clock-wise direction) represent light rain, shower, unorganized and organized storm events.

Figure 7 .
Figure 7. Overall structure of the sensitivity analysis test.

Figure 8 .
Figure 8. Influence of the a priori parameters on the efficiency of the algorithm in NSE criterion for light rain, shower, unorganized and organized storm events (clock-wise direction).The red dots in the figure indicate the location of the optimal parameters of both a priori and decorrelation distance.Note that, the NSE values are averaged over all maps in each event.

Figure 9 .
Figure 9. Example of rainfall maps retrieved in different events: (a) light rain, (b) shower, (c) organized and (d) unorganized storm.The figures on the left and right represent rainfall maps observed by weather radar (ground-truth) and retrieved by the proposed algorithm, respectively.

Figure 10 .
Figure 10.Comparison between inverse model and a priori knowledge at    ,    and 0.5.   resolutions which belong to 1st , 2nd and 3rd columns, respectively.Red solid lines indicate the NSE values using only a priori knowledge before applying the inverse algorithm.Blue solid lines indicate the NSE values obtained after applying the inverse algorithm.Rows (top 5

Figure 11 .
Figure 11.The retrieval efficiency in NSE at . .   resolution in different parts of the monitoring area classified as high, moderate, low and whole density.The evaluation test is in the presence of magnitude of measurement and model ( = 5 %) and 5

Figure 12 .
Figure 12.Scatter plot and fitted linear regression line (red) in 4 rain events: light rain, shower, organized and unorganized storm (top-down direction).The retrieved rainfall maps were obtained in the presence of the measurement and model (5 %), quantization (0.1 dB).Right and left hand side of the figure show the scatterplots computed for whole and high density regions of 5