Role of measurement uncertainty in the comparison of average areal rainfall methods

The main motivation for this research is the growing awareness of the impact of climate change and the increasing relevance of the United Nations Sustainable Development Goals, aiming to contribute to the measurement of quantities like precipitation and rate of rainfall. This knowledge is widely used in hydrology, climatology and meteorology, providing data and information applied in modelling, pattern definition and recognition, and forecasting. This work is concerned with estimating the average areal rainfall in a stipulated region from rainfall intensity observations made at measurement stations within that region. It focuses on three straightforward estimation approaches: the arithmetic mean method, the Thiessen polygon method and the isohyetal method. The evaluation of the associated measurement uncertainty, for which the law of propagation of uncertainty and a Monte Carlo method as described in guidance documents from the Joint Committee for Guides in Metrology are applied, is the main consideration. The approaches described may be readily applied by practitioners. A comparison of results from applying these methods to a simple example is made. Such results are required for conformity assessment and support in urban management and water resources management worldwide.


Introduction
The growing awareness of the impact of climate change and the United Nations Sustainable Development Goals [1] show the need to have reliable measurements of quantities to support urban and water resources management. Precipitation and rate From this definition (being the quantity interpreted as either the mass or volume of the liquid or solid products), precipitation intensity is a quantity defined as the amount of precipitation collected per unit time interval. The unit of precipitation is linear depth in mm (corresponding to a volume per unit area) and for liquid precipitation, kg m −2 (corresponding to a mass per unit area). The difference between rainfall and precipitation is that rainfall relates to water in its liquid state in the form of precipitated condensed droplets from atmospheric water vapour, while precipitation relates to the product of the condensation of atmospheric water vapour that falls under gravity. The measurement unit of rainfall intensity is linear depth per hour (mm h −1 ), usually obtained at a time interval of 1 min, or less in the case of extreme events or systems with high variability or intensity.
Observations of these quantities are often used to provide models for geographical surfaces, thus requiring rainfall and precipitation to be observed at various geographical locations. Difficulties in their measurement often arise from the influence of harsh conditions such as exposure, wind and topography, which have high impact on the accuracy of measurements. Wind effect is critical for the performance of instrumentation, leading to different shapes of gauges as shown in figure 1 [2], illustrating how streamlines of wind deformation are expected to affect the trajectory of precipitation particles, promoting a relevant error contribution to the measurement.
The recommended use of recording precipitation gauges is to have sufficient information related to the time scale and resolution to cater for the high variability of precipitation intensity. Such information is used in the technical process of reducing evaporation and wetting losses as sources of error. Three types of automatic precipitation recorders for measuring rainfall [3] are commonly used: 4 weighing-recording type, tipping-bucket type and floating type. The study carried out considered the use of the weighing and tipping-bucket types, both of which collect precipitation using an orifice and a funnel directed into receptacles allowing the volume to be weighed, in the first case, or multiples of volumes collected in a pair of buckets each having a reference volume quantity per second.
The WMO establishes reference conditions for the installation and use of these gauges [4], namely, the orifice height above the ground (commonly between 0.5 m and 1.5 m), conditions of the surroundings to avoid errors due to in-splashing from the ground, and specific geometries adopted for the orifice and the gauge. Wind field in the surroundings can be a major influence on the measurement. Special care is highly recommended by including windshields in the setup, establishing the type of surrounding surface, adopting the suggested relations with the vertical angle obstacles in the surroundings of the site, and choosing an appropriate gauge size and shape to minimize the wind effect [2].
To obtain the data needed for meteorology, climatology and hydrology predictive models, networks of stations are distributed in areas of interest in a way intended to represent the distribution of rainfall, being required to obtain rainfall intensity measurements at single points and combine them in order to calculate the volume of precipitation that falls over a given catchment area [5]. To achieve this aim, there are several methods that use sets of measurements at geographic locations and models to obtain the average areal rainfall.
In recent years, the increase in computational power and the use of geographic information systems [6] facilitated the development of complex hydrological models, resulting in improved measurement and forecast accuracy. Algorithmic implementations of these models, used to understand and to predict climatic phenomena, are highly dependent on extreme events, topography, accuracy of measurement instruments and quantity of climatic data [7]. It is necessary to find a balance between having a wide network of measurement stations and the related high costs of these network infrastructures.
The use of interpolation methods to estimate the spatial distribution of rainfall from rain gauge measurement data is a common technique that can be divided into two main groups [7]: A general approach for obtaining the volume of precipitation P g over a designated area given rainfall intensity values P i , i = 1, . . . , n, measured at n observation points, is as follows. Together with appropriate weights λ i , the P i are used to form the volume of precipitation from the linear combination According to a comparative analysis of techniques for spatial interpolation of precipitation [7], the methods most frequently used for this purpose are: • Arithmetic mean method, based on the average of observations for a set of stations, applied, for example, in flat areas; • Thiessen polygon method, also known as the nearest neighbour method [8], based on constructed polygon areas related to the nearest observation stations; • Inverse distance weighting method, based on weights defined as the inverse of the distance between observation stations, normalized appropriately; • Polynomial interpolation method, using an algebraic or a trigonometric polynomial function [9]; • Spline interpolation method, a mathematical model forming a minimum-curvature surface; • Moving window regression method, applying linear regression to the areas having a relationship between two variables (e.g. rainfall and elevation) [10].
These methods do not consider the contribution of topography, which affects the catchment of precipitation at a gauge. The isohyetal method addresses this contribution, allowing the drawing of isohyets (lines of equal rainfall depths designed to take into account the locations and the topography).
There are different types of kriging [7,11] considering how the mean value of the quantity of interest is applied: • Simple kriging, assuming the mean is constant and known; • Ordinary kriging, assuming the mean is constant but unknown; and • Universal kriging, assuming that the mean is represented by a polynomial function.
Some other types, for example, kriging for uncertain data [12,13], which take into account the variability due to rain gauge uncertainty, applying independently for the several observation points and introducing spatial and temporal variance into the modelling, are becoming more of interest. For this type of approach, recommended studies are for cases where rainfall observations have similar or different uncertainties and [14,15].
The use of kriging methods requires a spatially representative set of input data to achieve accurate results for areal rainfall. The study developed in the context of this work aims to use simple tools for practitioners to estimate areal rainfall and evaluated the associated measurement uncertainty in cases where there is a small group of observation stations, as is common practice in many regions. In such circumstances, users adopt relatively straightforward deterministic methods rather than those employing more complex data processing to produce their results.
For this study, the following spatial interpolation methods were considered: arithmetic mean method, Thiessen polygon method and isohyetal method. These methods give alternative interpretations of the physical quantity in the geometric context. For this reason, measurement uncertainty in the output plays a relevant role in the comparison of the accuracy of the methods. A brief description of each method for estimating the average areal rainfall is presented: (a) Arithmetic mean method: evaluates the arithmetic mean of considered single-point observations for a given area; (b) Thiessen polygon method: applies a graphical approach that defines relative polygonal areas related to each single station observations, providing a weighted sum of the observations; (c) Isohyetal method: applies a graphical approach based on the drawing of isohyet lines of equal rainfall, combining the observations weighted by the coverage areas between these lines.
As mentioned, the process used to obtain the average areal rainfall has two stages: the first obtains the rainfall intensity at each point of the network (gauge station); the second uses the information provided by the stations in the network to combine the point intensities using the above methods.

Rainfall intensity measurement model
The process initially involves the measurement of rainfall intensity, P i , at several locations, I, and the evaluation of the related measurement uncertainties. For this purpose, different types of techniques and rain gauges are used (tipping bucket, weighing, optic, etc), usually operating by collecting a quantity of liquid precipitation over a defined period [4]. For this study, two of the most common instruments were used, weighing gauges (figures 2 and 3) and tipping-bucket gauges (figures 4 and 5). In the first case, amounts of water related to weighing observations over time are collected, while, in the second case impulses are generated for a fixed volume.
Weighing gauges effectively are balances that record the mass of accumulated precipitation over time. The container should have a large capacity to account for the fact that this system does not empty itself. Solutions for its use in harsh climate conditions are also required (for example, the use of oil to avoid large evaporation effects, and antifreeze solutions).
The operation of tipping-bucket gauges allows water collection and guidance to a twin-bucket balance with both parts having an equal weight and reference volume. Every time a bucket   becomes full the balance changes position within a pivot axis, the other bucket moving into position to collect water while the first empties the collected water. The time between each change (tipping) is measured allowing the calculation of the rate of rainfall. Tipping-bucket gauges employ a contact closure (reed switch or relay contact), such that each tip produces an electrical impulse as a signal output. This output is recorded by a data logger or an ADC (analogue-to-digital converter, data acquisition system equipped with reed-switch reading ports). This mechanism provides continuous measurement without manual interaction.
The nature of the measurement of precipitation is affected by many natural conditions, implying the need to account for corrections and to evaluate the effect of errors in the methods [16,17]. Reports issued by WMO point out the need to use models to adjust the measured precipitation [18] based on corrections obtained from statistical data. Regarding errors (systematic effects), WMO also collected information provided by research, being able to state [4] 'the amount of precipitation measured by commonly used gauges may be less than the actual precipitation reaching the ground by up to 30% or more'. Considering the interest of this study in the rainfall intensity measurements obtained using tipping-bucket gauges and weighing gauges (other types like floating gauges and optical gauges were not considered for this purpose), data provided in the WMO reports [3,19] were taken into account. The assessment of errors in precipitation measurement usually relates its origin to the effects of wind, wetting and evaporation losses [18]. A general description of these sources is given in [4], including those based on [18]: (a) Error due to systematic wind field deformation above the gauge orifice: typically 2% to 10% for rain and 10% to 50% for snow; (b) Error due to the wetting loss on the internal walls of the collector; (c) Error due to the wetting loss in the container when it is emptied: typically 2% to 15% in summer and 1% to 8% in winter, for (b) and (c) together; (d) Error due to evaporation from the container (most important in hot climates): 0% to 4%; (e) Error due to blowing and drifting snow; (f ) Error due to the in-and out-splashing of water: 1% to 2%; (g) Systematic mechanical and sampling errors, and dynamic effects errors (i.e. systematic delay due to instrument response time): typically 5% to 15% for rainfall intensity, or even more in high-rate events (see 3); (h) Random observational and instrumental errors, including incorrect gauge reading times.
Considering these sources of error and uncertainty, a functional relation for precipitation (rain and snow contributions) was proposed by WMO [18], and adapted in 1990 by Legates and Willmott [20] as P k = k r P cr + k s P cs , where P cr = P gr + ΔP 1r + ΔP 2r + ΔP 3r + ΔP 4r , P cs = P gs + ΔP 1s The quantities in expression (1) are as follows: Subscript r-relates to liquid precipitation 'rain'; Subscript s-relates to 'solid' precipitation; P k -adjusted precipitation; k-adjustment factor for the effects of wind field deformation; P c -amount of precipitation caught by the gauge collector; P g -measured amount of precipitation in the gauge; ΔP 1 -adjustment for the wetting loss on the internal walls of the collector; ΔP 2 -adjustment for wetting loss in the container after emptying; ΔP 3 -adjustment for evaporation from the container; and ΔP 4 -adjustment for systematic mechanical errors. The adjustment factor k is a variable obtained from studies undertaken by Nespor and Sevruk [21], in which the ratio of correct to measured precipitation of rain and snow was studied using two unshielded gauges in different weather conditions of wind speed and intensity. The measurement of rainfall intensity, in units of mm h −1 , using weighing gauges or tippingbucket gauges, starts respectively with the measurement of mass or volume in units of time. The measurand is affected by sources of uncertainty according to the relational function (1), considering only the liquid contributions, P k = k r P cr = k r P gr + ΔP 1r + ΔP 2r + ΔP 3r + ΔP 4r . (2) The measurement uncertainty associated with an estimate of P k can be evaluated using the law of propagation approach described in the Guide to the expression of uncertainty in measurement (GUM) [22] or a Monte Carlo method (MCM) in GUM supplement 1 (GUM-S1) [23].

Average areal rainfall methods
The primary measurand, P av , is the average areal rainfall as a function of the rainfall intensity values or averages obtained at the m locations of measurement stations, being a function of a weighted linear combination of the P i : The function (3) is considered generic, being applied to all three methods studied.
In a general approach, the methods consider a number m of measurement stations providing rainfall values P i (i = 1, . . . , m) distributed across a basin. To illustrate the approach adopted in each method, as a starting point consider m = 4 measurements at rainfall stations located in a basin, as in figure 6.
The arithmetic-mean method evaluates the average of the estimates obtained at each location without establishing a relation between the position of the stations and the geometry of the area of observation. The evaluation considers that all weights are equal to 1/m: The second method is the method of Thiessen polygons, which uses a given set of locations to partition the plane into The average areal rainfall [equation (5)] using the Thiessen polygon method is the weighted mean the w i being the weights given by the relative areas of the polygons obtained, A i is the area of the polygon related to station i, and A is the total area of the basin. Any change of rainfall at the stations does not affect the geometry of the polygons. The third method studied uses weights proportional to contour map areas according to the location of isohyets (lines on a map or chart connecting areas of equal rainfall). The isohyetal method (figure 8) uses the single-point station information to establish contour map areas [24,25], with weights obtained by multiplying each contour area by the average rainfall in the area.
In the example (where m = 4), using the same basin and considering locations, m + 1 isohyets (P 1 , . . . , P m+1 ) are defined (P a , P b , P c , P d , P e in the figure), dividing the basin into m contour areas (A 1 , . . . , A m ), allowing the determination of the areal rainfall average, P av,isoh , using

Uncertainty propagation and average areal rainfall methods comparison
The propagation of uncertainty in this study of the average areal rainfall intensity has two stages. The first stage requires the evaluation of the uncertainty related to the measurement of rainfall intensity in the individual locations, using either  weighing gauges or tipping-bucket gauges. The second stage requires the evaluation of the uncertainty of average areal rainfall intensity using one of the three methods mentioned. Given the nature of precipitation phenomena and the variability inherent in the main sources of measurement error, the quantification of the resulting effects is usually difficult to establish. For the purpose of obtaining estimates of these quantities, references [26][27][28][29][30][31] were consulted. In the first stage, for both type of gauges, the input quantities are described in expression (2). A probability density function (PDF) was assigned to each quantity based on knowledge of that quantity. The mean of that PDF was used as an estimate of the quantity and the standard deviation as the associated standard uncertainty. The input quantities are shown in table 1 together with the assigned PDFs and their relative standard uncertainties (in the table the index 'r' is suppressed).
Since the PDFs for wetting loss are not symmetrical about zero, a correction was made to the estimate equivalent to the half width of the PDF and a zero-centred PDF was used in the uncertainty evaluation.
The second stage accounts for the uncertainty associated with the measurement of rainfall intensity at each location, for both types of gauges considered, having as input the combined average areal rainfall uncertainty obtained in stage 1. In this case, the evaluation of measurement uncertainty does not account for possible correlation between measurements at the different locations.
The evaluation of the uncertainty for a measurement of rainfall precipitation of 10 mm h −1 , was made using RStudio [32], with 10 7 Monte Carlo trials for each calculation. The evaluation allowed, for both gauges, the PDF for the output quantity, P k , and the relative expanded uncertainty, U 95 (P k ), for a confidence interval of 95%, by applying GUM [22] and GUM-S1 [23], to be provided. The values obtained and the related PDFs are shown in table 2 and in figures 9 and 10 for the weighing gauge and tipping-bucket gauge, respectively. These figures also show the scaled histograms produced using GUM-S1 and used as a basis for the (continuous) PDFs shown in blue.
The results show consistency with the normal distribution in the case of the weighing gauge and a small deviation from normality in the case of tipping-bucket gauge, identified by skewness and kurtosis parameter values that differ slightly from normal reference values of 0 and 3, respectively. Table 1. Input quantities, relative standard uncertainties and assigned PDFs related to rainfall intensity measurement using weighing gauges and tipping-bucket gauges.

Quantity Description
Weighing gauge standard uncertainty/% Tipping-bucket gauge standard uncertainty/% PDF k Error due to systematic wind field deformation above the gauge orifice 5 5 Gaussian P g1 Error due to the in-and out-splashing of water 2 2 Uniform P g2 Random observational and instrumental errors 2 2 Uniform ΔP 1 Error due to the wetting loss on the internal walls of the collector n/a 5 Uniform a ΔP 2 Error due to wetting loss in emptying the container 1 5 Uniform a ΔP 3 Error due to evaporation from the container 2 1 Uniform ΔP 4 Systematic mechanical and sampling errors, and dynamic effects errors 5 10 Uniform a These uniform PDFs have intervals from 0% to 5% (asymmetric with respect to 0%), considering that loss quantity would increase the estimate and negative values are not achievable. Table 2. Parameters and expanded measurement uncertainties obtained for the weighing gauge and tipping-bucket gauge using GUM and GUM-S1 (CI 95 = 95% coverage interval).
Weighing gauge /mm h −1 Tipping-bucket gauge /mm h −1 GUM GUM-S1 GUM GUM-S1 Skewness (GUM-S1) 0.08 Skewness (GUM-S1) 0.11 Kurtosis (GUM-S1) 2.96 Kurtosis (GUM-S1) 2.73 For the second stage, a comparison of the measurement uncertainty associated with the average areal rainfall is made for the three methods considering that local measurement was made using either weighing or tipping-bucket gauges.
For this study, figure 6 was adopted as being representative of a territory having four measurement stations with the following daily average rainfall: P 1 = 12 mm, P 2 = 18 mm, P 3 = 36 mm and P 4 = 28 mm. The relative standard measurement uncertainty considered for each estimate was based on the evaluation obtained at stage 1: 6% for weighing gauges and 8% for tipping-bucket gauges.
The first approach to calculate the daily average areal rainfall used the arithmetic mean method, applying equation (4). To evaluate the measurement uncertainty, u (P av ), for this linear model, the GUM approach gives an exact solution, assuming that measurements at the various stations are uncorrelated: Figure 10. Tipping-bucket gauge PDFs obtained using GUM (red line) and GUM-S1 (blue line) and scaled histogram of output numerical sequence.
In the given catchment area, using the values given above for P 1 to P 4 , the estimate of P av (daily average areal rainfall) is 23.5 mm. Considering that u (P i ) = 0.06P i for weighing gauges and u (P i ) = 0.08P i for tipping-bucket gauges, expression (7) is used to obtain the standard uncertainties for both types of gauges considered in stage 1: u (P av ) = 0.76 mm, using weighing gauge standard uncertainty, u (P av ) = 1.0 mm, using tipping-bucket gauge standard uncertainty.
The second method, the Thiessen polygon method, was applied using the same P i values but it requires the evaluation of the areas of the polygons that gives the weights considered in expression (5). The areas related to the polygons shown in figure 9 were obtained using a planimeter technique, giving approximate values of To evaluate measurement uncertainty using the Thiessen polygon method, uncertainty contributions for daily average areal rainfall were taken to be the same as for the arithmetic mean method. The combined uncertainty also takes account of the uncertainty contributions related to the area weight of each polygon, u(A i /A), assessed to be 0.01. In this case, the evaluation of the measurement uncertainty associated with the daily average areal rainfall intensity used a MC approach according to GUM-S1. The numerical evaluation was developed for both types of gauges, using RStudio, with 10 6 MC trials. Using expression (5) the daily average areal rainfall, P av,Tp , is 21.3 mm and u P av·Tp ≈ 0.8 mm, using weighing gauge standard uncertainty, u P av·Tp ≈ 1.0 mm, using tipping-bucket gauge standard uncertainty.
The third approach applies the isohyetal method to the same area, requiring the values for the isohyets presented in figure 8, in order to make the computation of the annual average areal rainfall intensity according to expression (6). In this case, the values considered for the isohyets, considering the average estimates of P 1 to P 4 , were the following: P a = 6 mm; P b = 15 mm; P c = 24 mm; P d = 32 mm; P e = 44 mm.
Standard measurement uncertainties considered for the isohyets were taken as those used previously. As in the second method, the relative areas between adjacent isohyets need to be evaluated, obtaining The standard uncertainty related to the area weight of each sub-area, u(A i /A), were taken to be 0.01.
The evaluation of the measurement uncertainty associated with the daily average areal rainfall intensity again used a MC approach according to GUM-S1. The numerical evaluation was developed for both types of gauges, using RStudio, with 10 7 runs for each calculation. Using expression (6) the estimate of P av,Isoh is 21.9 mm and the associated standard uncertainties are u (P av·Isoh ) ≈ 1.9 mm, using weighing gauge uncertainty, u (P av·Isoh ) ≈ 2.4 mm, using tipping-bucket gauge uncertainty.
A summary of the results is given in tables 3 and 4 and the output PDFs in figures 11 and 12, respectively, for the use of weighing gauges and tipping-bucket gauges as measurement instruments.  Figure 11. Comparison of results and 95% expanded uncertainty, using as input the uncertainty rainfall intensity measurement of weighing gauges, for the arithmetic mean method (blue line), Thiessen polygon method (red line) and isohyetal method (black line).

Interpretation of comparison results
The diversity of measuring instruments has impact on the estimates and associated uncertainties, and the decisions taken should consider the effects due to the uncertainty contributions. In this study, two techniques (weighing and tippingbucket) for the same measurements (precipitation and rainfall) were adopted. In both cases the results obtained using GUM and MCM (GUM-S1) were consistent. In the case of the tipping-bucket gauge, a slightly higher degree of flatness of the PDF for rainfall precipitation was detected using GUM-S1, indicating a small deviation from normality. In both cases, average estimates were corrected to consider the contribution of the wetting loss bias, which could not be measured in the process.
Regarding the three methods considered, commonly applied to evaluate daily, monthly or annual average areal rainfall, the interest of the studies carried out were related to performing a comparison of the results based on the measurement uncertainties associated with the estimates provided. The results obtained and presented in tables 3 and 4 and in figures 11 and 12 show differences in the estimates of the average areal rainfall, from 21.3 mm to 23.5 mm. Regarding the 95% expanded measurement uncertainty obtained (tables 3 and 4), agreement was found between the arithmetic mean method and Thiessen polygon method with the measurement uncertainty around 20% higher for the isohyetal method. The comparison between weighing gauges and tipping-bucket gauges showed a difference of 30%, enhancing the conclusions that the impact of the type of gauge used and the method adopted are high.
The studies carried out were able to show some interesting features of the models and the way they affect the measurand of interest. Further studies should be made to include the effect of correlation that was not considered in this simple analysis, and the effect of conditions related with the dynamics of the measurement process and the use of corrective algorithms related to post-processing of data [32][33][34][35].