Ensemble-based statistical interpolation with Gaussian anamorphosis for the spatial analysis of precipitation

Hourly precipitation over a region is often simultaneously simulated by numerical models and observed by multiple data sources. An accurate precipitation representation based on all available information is a valuable result for numerous applications and a critical aspect of climate. Inverse problem theory offers an ideal framework for the combination of observations with a numerical model background. In particular, we have considered a modified ensemble optimal interpolation scheme, that takes into account deficiencies of the background. An additional source of uncertainty for the ensemble background has 5 been included. A data transformation based on Gaussian anamorphosis has been used to optimally exploit the potential of the spatial analysis, given that precipitation is approximated with a gamma distribution and the spatial analysis requires normally distributed variables. For each point, the spatial analysis returns the shape and rate parameters of its gamma distribution. The Ensemble-based Statistical Interpolation scheme with Gaussian AnamorPhosis (EnSI-GAP) is implemented in a way that the covariance matrices are locally stationary and the background error covariance matrix undergoes a localization process. Con10 cepts and methods that are usually found in data assimilation are here applied to spatial analysis, where they have been adapted in an original way to represent precipitation at finer spatial scales than those resolved by the background, at least where the observational network is dense enough. The EnSI-GAP setup requires the specification of a restricted number of parameters and specifically the explicit values of the error variances are not needed, since they are inferred from the available data. The examples of applications presented provide a better understanding of the characteristics of EnSI-GAP. The data sources considered 15 are those typically used at national meteorological services, such as local area models, weather radars and in-situ observations. For this last data source, measurements from both traditional and opportunistic sensors have been considered. Copyright statement. Usage rights are regulated through the Creative Commons Attribution 3.0 License (https://creativecommons.org/ licenses/by/3.0).

covariance matrix is derived from numerical model ensemble and where Gaussian anamorphosis is applied directly to precipitation data. An additional innovative part of the method is that EnSI-GAP does not require the explicit specification of error 95 variances for the background or observations, as most of the other methods (Soci et al., 2016). In fact, those error variances are often difficult to estimate in a way that is general enough to cover a wide range of cases. Our approach is to specify the reliability of the background with respect to observations, in such a way that error variances can vary both in time and space. An additional innovative part of our research is that we consider opportunistic sensing networks of the type described by de Vos et al. (2020) within the examples of applications proposed. Thanks to those networks, for some regions we can rely on an 100 extremely dense spatial distribution of in-situ observations. The remaining of the paper is organized as follows. Sec. 2 describes the EnSI-GAP method in a general way, without references to specific data sources. Sec. 3 presents the results of EnSI-GAP applied to three different problems: an idealized situation, then two examples where the method is applied to real data, such as those mentioned above. The results are discussed in Sec. 3.3.

2 Methods: EnSI-GAP, Ensemble-based statistical interpolation with Gaussian anamorphosis for precipitation
We assume that the marginal probability density function (PDF) for the hourly precipitation at a point in time follows a gamma distribution (Wilks, 2019). This marginal PDF is characterized through the estimation of the gamma shape and rate for each point and hour.
Precipitation fields are regarded as realizations of locally-stationary, trans-Gaussian random fields, where each hour is con-110 sidered independently from the others. Trans-Gaussian random fields are used for the production of precipitation observational gridded datasets by Frei and Isotta (2019). A random field is said to be stationary if the covariance between a pair of points depends only on how far apart they are located from each other. Precipitation totals are nonstationary random fields because the covariance between a pair of points in space depends not only on the distance between them but it varies also when considered in different directions. In our method, precipitation is locally modeled as a stationary random field. The covariance parame-115 ter estimation and spatial analysis are carried out in a moving-window fashion around each grid point. A similar approach is described by Kuusela and Stein (2018) and the elaboration over the grid can be carried out in parallel for several grid points simultaneously.
A particular implementation of EnSI-GAP is reported in Algorithm 1. The mathematical notation and the symbols used are described in two tables: Tab. 1 for global variables and Tab. 2 for local variables, which are those variables that vary from 120 point to point. As in the paper by Sakov and Bertino (2011), upper accents have been used to denote local variables. If X is a matrix, X i is its ith column (column vector) and X i,: is its ith row (row vector). The Bayesian statistical method used in our spatial analysis is optimal for Gaussian random fields. Then, a data transformation is applied as a pre-processing step before the spatial analysis. The introduction of a data transformation compels us to inverse transform the predictions of the spatial analysis back into the original space of precipitation values. Algorithm 1 can be divided into three parts, that are described in 125 the next sections: the data transformation in Sec. 2.1, the Bayesian spatial analysis in Sec. 2.2 and the inverse transformation in Sec. 2.3.

Data transformation via Gaussian anamorphosis
The data transformation chosen is a Gaussian anamorphosis (Bertino et al., 2003), that transforms a random variable following a gamma distribution into a standard Gaussian. This pre-processing strategy has been used in several studies in the past, e.g. Amezcua and Leeuwen (2014); Lien et al. (2013). A visual representation of the transformation process can be found in Fig.1 of the paper by Lien et al. (2013).
The hourly precipitation background and observations,X f andỹ o respectively, are transformed into those used in the spatial analysis by means of the Gaussian anamorphosis g(): The gamma shape and rate are derived from the precipitation values by a fitting procedure based on maximum likelihood.
The maximum likelihood estimators are calculated iteratively by means of a Newton-Raphson method as described by Wilks (2019), Sec. 4.6.2. In particular, the gamma distribution parameters are fitted to each ensemble member field of precipitation separately. Then, the averaged shape and rate are used in g() for Eqs. (1)-(2).

140
In Gaussian anamorphosis, zero precipitation values must be treated as special cases, as explained by Lien et al. (2013). The solution we adopted is to add a very small amount to zero precipitation values, e.g. 0.0001 mm, then to apply the transformation g() to all values. The same small amount is then subtracted after the inverse transformation. This is a simple but effective solution for spatial analysis, as shown in the example of Sec. 3.1.1. In principle, the statistical interpolation is sensitive to the small amount chosen, such that using 0.01 mm instead of 0.0001 mm will return slightly different analysis values in the 145 transition between precipitation and no-precipitation. In practice, we have tested it and we found negligible differences when values smaller than e.g. 0.05 mm (half of the precision of a standard rain gauge measurement) have been used.

Spatial analysis
The spatial analysis inside Algorithm 1 has been divided into three parts. In Sec. 2.2.1, global variables have been defined. Then, as stated in the introduction of Sec. 2, the analysis procedure is performed on a gridpoint-by-gridpoint basis. In Sections 2.2.2-150 2.2.3, the procedure applied at the generic ith gridpoint is described. In Sec. 2.2.2, the specification of the local error covariance matrices is described. In Sec. 2.2.3, the standard analysis procedure is presented together with the treatment of a special case.

Definitions
In Bayesian statistics, according to Savage (1972), a state is "a description of the world, which is the object which we are concerned, leaving no relevant aspect undescribed" and "the true state is the state that does in fact obtain". The object of our 155 study is the hourly precipitation field x(), that is the hourly total precipitation amount over a continuous surface covering a 5 https://doi.org/10.5194/npg-2020-20 Preprint. Discussion started: 19 June 2020 c Author(s) 2020. CC BY 4.0 License. spatial domain in terrain-following coordinates r. Our state is the discretization over a regular grid of this continuous field. The true state (our "truth", x t ) at the ith grid point is the areal average: where V i is a region surrounding the ith grid point. The size of V i determines the effective resolution of x t at the ith grid point.

160
Our aim is to represent the truth with the smallest possible V i . The effective resolution of the truth will inevitably vary across the domain. In observation-void regions, the effective resolution will be the same as that of the numerical model used as the background, then approximately o(10 − 100 km 2 ) for high-resolution local area models (Müller et al., 2017). In observationdense regions, the effective resolution should be comparable to the average distance between observation locations, with the model resolution as the upper bound.

165
The analysis is the best estimate of the truth, in the sense that it is the linear, unbiased estimator with the minimum error variance. The analysis is defined as x a = x t + η a , where the column vector of the analysis error at grid points is a random variable following a multivariate normal distribution η a ∼ N (0, P a ). The marginal distribution of the analysis at the ith grid point is a normal random variable and our statistical interpolation scheme returns its mean value x a i and its standard deviation σ a i = P a ii .

170
As for linear filtering theory (Jazwinski, 2007), the analysis is obtained as a linear combination of the background (a priori information) and the observations. The background is written as x b = x t +η b , where the background error is a random variable The background PDF is determined mostly, but not exclusively, by the forecast ensemble, as described in Sec. 2.2.1. The forecast ensemble mean is x f = k −1 X f 1, where 1 is the m-vector with all elements equal to 1. The background expected value is set to the forecast ensemble mean, plays a role in the determination of P b , as defined in Sec. 2.2.2.
The p observations are written as y o = Hx t + ε o , where the observation error is ε o ∼ N (0, R) and H is the observation operator, that we consider as a linear function mapping R m onto R p . 180

Specification of the observation and background error covariance matrices
Our definitions of the error covariance matrices follow from a few working assumptions, WAn indicates the nth working assumption and the abbreviations will be used in the text. WA1, background and observation uncertainties are weather-and location-dependent. WA2, the background is more uncertain where either the forecast is more uncertain or observations and forecasts disagree the most. WA3, observations are a more accurate estimate of the true state than the background. We 185 want to specify how much more we trust the observations than the background in a simple way, such as e.g. "we trust the observations twice as much as the background". WA4, the local observation density must be used optimally to ensure a higher effective resolution, as it has been defined in Sec. 2.2.1, where more observations are available. WA5, the spatial analysis at a particular hour does not require the explicit knowledge of observations and forecasts at any other hour. However, constants in the covariance matrices can be set depending on the history of deviations between observations and forecasts. WA5 makes the 190 procedure more robust and easier to implement in real-time operational applications.
A distinctive feature of our spatial analysis method is that the background error covariance matrix i P b is specified as the sum of two parts: a dynamical component and a static component. The dynamical part introduces nonstationarity, while the static part describes covariance stationary random variables. This choice follows from WA1 and it has been inspired by hybrid data assimilation methods (Carrassi et al., 2018). The dynamical component of the background error covariance matrix is obtained 195 from the forecast ensemble. Because the ensemble has a limited size, and often the number of members is quite small (order of tenths of member), a straightforward calculation of the background covariance matrix will include spurious correlations between distant points. Localization is a technique applied in DA to fix this issue (Greybush et al., 2011). The static component has also been introduced to remedy the shortcomings of using numerical weather prediction as the background. There are deviations between observations and forecasts that cannot be explained by the forecast ensemble. A typical example is when 200 all the ensemble members predict no precipitation but rainfall is observed. In those cases, we trust observations, as stated through WA3. Then, the static component adds noise to the model-derived background error, as in the paper by Raanes et al. (2015). In Bocquet et al. (2015), the static component is referred to as a scale matrix, since it is used to scale the noise component of the model error, and we adopt the same term here. We will also refer to this matrix, and its related quantities, with the letter u to emphasize that this component of the background error is "unexplained" by the forecast.

205
i P b is written as: The first component on the right-hand side of Eq. (5) is the dynamical part. P f is the forecast uncertainty of Eq. (4), i Γ is the localization matrix and • is the Schur product symbol. The localization technique we apply is a combination of local analysis and covariance localization, as they have been defined by Sakov and Bertino (2011). In the local analysis, only the functions have been described for instance by Gaspari and Cohn (1999). We have chosen not to inflate or deflate P f directly and to modulate the amplitude of background covariances only through the terms of Eq. (5), this way we reduce the number of parameters that need to be specified. As a matter of fact, for the combination of observations and background in the analysis 220 procedure, the m by m covariance matrices are never directly used. Instead, the matrices used are: the covariances between f ≥ i σ 2 ob /(1 + ε 2 ) and when the forecast ensemble spread is too small, then Eq. (11) ensures that Eq. (10) is valid by setting

Analysis procedure
The expressions for the analysis and its error variance are direct results of the linear filter theory (Jazwinski, 2007)  variance is derived from the error covariance matrices: Eq. (13)-(14) are also typical of Optimal Interpolation and the formulation used is similar to the one adopted by Uboldi et al.

275
(2008), which follows from Ide et al. (1997). Because an ensemble of fields is used to determine the background error, the method is similar to Ensemble OI (Evensen, 1994).
The special case of i σ 2 f = 0 and i σ 2 ob /(1 + ε 2 ) = 0 still needs to be considered, since it does not belong to the cases of Eqs. (11)-(12). This is the case when both observations and forecast ensemble means at observation locations have exactly the same values and, in addition, all the ensemble members have the same values too. In practice, this might happen only when no 280 precipitation is observed and forecasted at the same time. Since the innovation is a vector of zeros, from Eq. (13) it follows that Because all the information available shows an exceptional level of agreement, we have chosen to set the analysis error variance to zero, such that for those points the analysis PDFs are Dirac's delta functions.

Data inverse transformation
The direct inverse transformation at the ith grid point is written as: however, we need to back-transform a Gaussian PDF and not a scalar value. Eq. (15) returns an estimate of the median of the gamma distribution associated to the ith grid point. Our goal is to obtain the gamma shape and rate. To achieve that, the direct inverse transformation g −1 is applied to 400 quantiles of the univariate gaussian PDF defined byx a i and (σ 2 ) a i , a similar approach is used by Erdin et al. (2012). Then, a least-mean-square optimization procedure is used to obtain the optimal shape 290 and rate that better fits the back-transformed quantiles. Given the optimal shape and rate, it is possible to obtain the statistics that better represent the distribution for a specific application. In Sec. 3, the analysis value chosen is the mean as it is the value that minimizes the spread of the variance. However, other choices may be more convenient depending on the applications, as discussed by Fletcher and Zupanski (2006) where, for instance, the mode was chosen as the best estimate. In Sec. 3, we will also consider selected quantiles of the gamma distribution to represent analysis uncertainty.

Results
EnSI-GAP is applied to two case studies in Sec. 3.1, one on synthetic and one on real-world data. In addition, its performances are evaluated over several months in Sec. 3.2. Then, a discussion of the results is presented in Sec. 3.3. In Sec. 3.1.1, EnSI-GAP is applied over a one-dimensional grid and in a "controlled environment", that is over synthetic data specifically generated for testing EnSI-GAP on precipitation. The spatial analysis is performed with and without data 300 transformation to assess its effects. Furthermore, we have compared different ways of specifying the scale matrix and we have investigated the sensitivity of EnSI-GAP to variations of α and ε 2 . In Sec. 3.1.2, a second more realistic example of application for EnSI-GAP is reported, where the spatial analysis is performed for a case study of intense precipitation, which happened in July 2019 over South Norway. The data are those used in the operational daily routine at MET Norway. The forecasts are from the MetCoOp Ensemble Prediction System (MEPS, 305 Frogner et al., 2019). MEPS has been running operationally four times a day (00 UTC, 06 UTC, 12 UTC, 18 UTC) since November 2016 and its ensemble consists of 10 members. The hourly precipitation fields are available over a regular grid of 2.5 km. The observational dataset of hourly precipitation is composed of two data sources: precipitation estimates derived from the composite of MET Norway's weather radar and meteorological weather stations equipped with ombrometers, such as rain gauges or other devices. The hourly precipitation in-situ observations have been retrieved from MET Norway's climate 310 database frost.met.no (last time the website was checked is 2020-05-13). In addition to MET Norway's official weather stations, the database includes data collected by several Norwegian public institutions, such as for example: universities (e.g., the Norwegian Institute of Bioeconomy Research -Nibio), the Norwegian Water Resources and Energy Directorate (NVE), the Norwegian Public Roads Administration (Statens vegvesen). As described in the recent paper by Nipen et al. (2020), MET Norway is successfully integrating amateur weather stations temperature data into its operational routine. In this study, hourly The domain, data sources and grid settings of the cross-validation experiments are the same as for the case study of intense precipitation of Sec. 3.1.2. We are aiming at evaluating EnSI-GAP from both a deterministic and probabilistic viewpoint.
The verification scores considered are commonly used in forecast verification and described in several books, such as for example Jolliffe and Stephenson (2012). A further useful reference for the scores is the website of the World Meteorological 325 Organization https://www.wmo.int/pages/prog/arep/wwrp/new/jwgfvr.html (last time checked 2020-05-13).

One-dimensional simulations
EnSI-GAP is shown over a one-dimensional grid with 400 points and spacing of 1 spatial unit, or 1 u, such that the domain covers the region from 0.5 u to 400.5 u and the generic ith grid points is placed at the coordinate i u. A simulation begins 330 with the creation of a true state, then observations and ensemble background are derived from it. The statistical interpolation scheme is applied with different configurations in order to illustrate the behaviour of the method.
The simulation presented here is shown in Fig. 1a. For each grid point, the true value (black line) is generated by a random extraction from the gamma distribution with shape and rate set to 0.2 and 0.1, respectively. To ensure spatial continuity of the truth, an anamorphism is used to link a 400-dimensional multivariate normal (MVN) vector with the gamma distribution. The 335 samples from the MVN distribution, with a prescribed continuous spatial structure, are obtained as described by Wilks (2019), chapter 12.4. The MVN mean is a vector with 400 components all set to zero and the covariance matrix is determined using a Gaussian covariance function with 10 u as the reference length used for scaling distances. The effective resolution (Sec. 2.2.1) of the truth is then 10 u.
The ensemble background (gray lines in Fig. 1a) on the grid is obtained by perturbing the truth. The background values at 340 observation locations are obtained using nearest neighbour interpolations applied to the ensemble members. The observation operator H (Sec. 2.2.1) is the nearest neighbour interpolation, that is a matrix with all zeros except for a single element on each line equal to 1 at the element corresponding to the closest grid point. The truth is perturbed considering four main error sources, which are typically found in precipitation fields simulated by numerical models. The first typical error is the misplacement of precipitation events, that is implemented here independently for each ensemble member by shifting the true values along the 345 grid by a random number between -10 u and +10 u. Second, the effective resolutions of the background members are set to be coarser than the truth. For each member, the coarser resolution is obtained by multiplying the true values by a coefficient derived from a uniform distribution with values between 0.05 and 2 and a spatial structure function given by a MVN with Gaussian covariance function with a reference length greater than that of the truth, which is 10 u. The exact reference length scale varies between each member as it is extracted from a Gaussian distribution with mean of 50 u and a standard deviation 350 of 5 u. Third, as previously stated in Sec. 2, the challenging special case of the background showing no-precipitation while the true state reports precipitation is considered. For this purpose, all the ensemble members for the grid points between 200 u and 300 u are set to 0 mm. The fourth typical error is a variation of the third one, the background at grid points between 50 u and 150 u follows an alternative truth but different from 0 mm. It is possible to recognize the four regions described above on the grid in Fig. 1 by means of the coordinate on the x-axis. Because we had to ensure continuity of the background, we have 355 enforced smooth transitions between the regions mentioned above and their surroundings.
The number of observations (blue dots in Fig. 1a) is set to 40. The observation locations are randomly chosen, such that the observation density is variable across the grid. The observed value at a location is obtained as the true value of the nearest grid point, plus a random noise that is determined as a random number between -0.02 and 0.02 that multiplies the true value.
The procedure is consistent with the fact that observation precipitation errors should follow a multiplicative model (Tian et al.,360 2013). The observation distribution is denser in the central part of the domain and sparser closer to the borders. The constraints on the number of observations are: 5 between 1 u and 100 u; 30 between 101 u and 300 u; 5 between the 301 u and 400 u.
The panel d in Fig. 1 shows the Integral Data Influence (IDI Uboldi et al., 2008), that is a parameter that stay close to 1 for observation-dense regions, while it is exactly equal to 0 in observation-void regions. In practice, the IDI at the ith grid point is computed here as the analysis in Eq. (13) when all the observations are set to 1 and the background to 0, moreover only the Where the IDI is close to zero, the analysis is as good as the background. The simulation has been configured such that the domain is well covered by the observations. In this way, it can provide insights on the combination process itself. As it follows from Eq. (13) and Eq. (5), IDI depends on the error covariance matrices. We have used the equations of Algorithm 1, then the only additional parameter we have to set is D i . In this example, D i is estimated as the distance between the ith grid point and its closest 3rd observation location. This way, IDI stays close to 370 1, meaning that the observations do have an impact on the analysis, even for observation-sparse regions. D i is shown in Fig. 1c and it has been constrained to vary between 5 u and 20 u. More details on how to set D i are discussed in Sec. 3.3.
The effect of the Gaussian anamorphism is shown in Fig. 1b. The transformed precipitation varies within a smaller range than the original precipitation, thus effectively shortening the tail of the distribution, reducing its skewness and making it more similar to a Gaussian distribution.

375
The transformed precipitation analysis is shown in Fig. 2 for different configurations of the statistical interpolation scheme.
The localization matrix i Γ, of Eq. (5), for all panels is specified using Gaussian functions, of the form of those used in Algorithm 1 for i Z and i V, with L i constant and set to 25 u for all the grid points. The value of ε 2 is set to 0.1, which means that we trust much more the observations than the background. The parameters that are allowed to vary are those determining the scale matrix. As stated above, D i is determined adaptively at each grid point as shown in Fig. 1c. α is set to two different values. In 380 12 https://doi.org/10.5194/npg-2020-20 Preprint. Discussion started: 19 June 2020 c Author(s) 2020. CC BY 4.0 License. the left column (panels a and c) α = 0.1, while in the right column (panels b and d) α = 1. The background uncertainty is a trade off between the ensemble spread and the averaged innovation, with α = 1 the weight of each component is determined by the simulation presented in Fig. 1, without considering that is just one possible realization of the ensemble that should be used to derive robust parameter estimates. In the case of α = 0.1, the scale matrix is multiplied by a smaller value of i σ 2 u , as can be seen from Eq. (11). The panels in the top row (panels a and b) are obtained using a Gaussian function for i Γ u , as in Algorithm 1, 385 while the bottom row (panels c and d) shows the analyses when an exponential function is used. The exponential correlation function for the spatial analysis of precipitation is used for example by Mahfouf et al. (2007); Lespinas et al. (2015); Soci et al. (2016).
In Fig. 2, the analysis mean (or expected value, red line) fits the observations and the observed value is often within the analysis PDF envelope shown (pink region), despite the background deficiencies. The same comments hold true also when 390 ε 2 = 1 (not shown here), with an increase of the analysis residuals (analysis minus observation) and a slightly larger spread of the analysis PDF. The analysis means, in all panels, follow the observed values closely and variations in the scale matrix seem to have limited effects on the analysis means, at least in regions with a dense observational network, such as between 100 u and 300 u. The analysis ensemble spread is more sensitive to variations in the scale matrix. In the case of α = 1, the spread is larger than for α = 0.1 because of the increased background uncertainty. In all panels, the analysis presents the higher 395 uncertainties between 50 u and 100 u, where observations are sparse and the background follows an alternative truth. The analysis uncertainty is large also between 200 u and 300 u, where the background ensemble is set to no-precipitation while the rather dense observational network reports yes-precipitation. Interestingly enough, the region between 200 u and 300 u is also where the analyses derived from Gaussian and exponential functions differ the most. In this region, the ensemble-dependent with the data transformation, the spread associated to the analysis PDF is proportional to the observed value, moreover the skewness of the PDF allows for values higher than the mean to happen with a higher probability than in the case of no data transformation. Those two characteristics are positive achievements of EnSI-GAP, since they are typical of a gamma PDF that is regarded as a good statistical model for precipitation uncertainties (Wilks, 2019). The impact of the data transformation seems to be less significant over the analysis means, which, given the EnSI-GAP settings used, are principally constrained to fit 420 the observations. The panels a-d of Fig. 3 illustrate how the analysis uncertainty is modified when passing from the transformed to the original space. The regions with the higher uncertainties in Figs. 2-3 coincide, except for the noticeable difference of the points between 250 u and 300 u. In Fig. 2, the analysis PDFs for all those grid points show Gaussian functions with some spread. However, the same points in Fig. 3 show a much narrower analysis PDF, often without any significant spread even though we are not in the special case of For all panels, α = 1, L = 25 u and Gaussian covariance functions are used for localization and the specification of the scale matrix. D i is the same shown in Fig. 1c. The effect on the error variances of using a data transformation against not using it is investigated here. EnSI-GAP variances are shown in panels a and b, with ε 2 = 0.1 and ε 2 = 1 respectively. When ε 2 = 1, i σ 2 b is 435 equal to i σ 2 o everywhere but their values vary along the one dimensional grid. The scheme without data transformation is shown in panels c and d for the same ε 2 values. As a consequence of the data transformation, the scales on the y-axis differ between panels a-b and c-d. For all panels, the analysis PDF spread is higher between 50 u and 100 u, as can be seen e.g. in Fig. 3b that correspond to Fig. 4a. When the data transformation is applied, i σ 2 b reaches its maximum between 200 u and 300 u, with a second maximum between 50 u and 105 u. In the case without data transformation, the two maxima of i σ 2 b are still in the 440 same regions, though the principal one is between 50 u and 150 u. Besides, when using a data transformation the two maxima do have rather similar values, while without using data transformations the highest maximum can be four times the smallest one (Fig. 4d). In all panels, i σ 2 u is different from zero in the region between 180 u and 300 u. In addition, for ε 2 = 0.1, i σ 2 u is different from zero also in the region between 50 u and 150 u because the variation of ε 2 = 0.1 modifies the threshold value of i σ 2 ob /(1 + ε 2 ) in Eqs. (11)-(12). Then, one interpretation is that between 200 u and 300 u, when the background fails to predict 445 the occurrence of precipitation, all spatial analysis scheme recognize the importance of using the additional term modulated by i σ 2 u for the specification of the background error. On the other hand, for the type of background failure taking place between 50 u and 150 u, that is the background is following an alternative truth with respect to observations, the behaviour of the spatial analysis scheme is less predictable and it is more sensitive to the EnSI-GAP settings. Furthermore, this example indicates that a data transformation is useful to keep the variances at reasonable values.

Intense precipitation case over South Norway
A mass of moist air from the ocean moving towards the Norwegian mountains was at the origin of several intense showers over western Norway on the 30th of July 2019. South Norway is the domain under consideration and it is shown in Fig. 5, the domain measures 373 km in the meridional direction and 500 km in the zonal direction. The measurements from the weather stations managed by MET Norway show hourly precipitation values with more than 20 mm, which is extremely intense given 455 the climatology of the region. In addition, thousands of lightning strikes have been recorded (not shown here), thus confirming the convective nature of the precipitation. Intense events have been observed in the afternoon along the coast and over the nearby mountains, especially in Sogn og Fjordane. This region is shown as the black box in Fig. 5, it extends for 80 km in both meridional and zonal directions. Point A corresponds to the center of a grid box where a maximum of precipitation has been observed. Point B is the center of a grid box that is not covered by observations and where a maximum of precipitation has 460 been reconstructed by the analysis. The distance between points A and B is 14 km. In Sogn og Fjordane damages have been reported, they were caused by the heavy rain that also triggered a series of landslides. One of them caused a fatality when a driver was caught in the debris flow.
On both domains, the focus is on the representation of hourly precipitation patterns at the mesoscale, as it has been defined by Thunis and Bornstein (1996); Stull (1988), though over different domains we will focus on different parts of the mesoscale.

465
South Norway is used to show that the variability of the fields represented by the forecast ensemble members involves mostly the Meso-β part of the mesoscale (i.e. spatial scales from 20 km to 200 km). Sogn og Fjordane is a domain where highresolution information is needed to support fine-scale analysis by e.g. civil protection authorities. In this case, we will study precipitation patterns at the Meso-γ scale (i.e. from 2 km to 20 km). The EnSI-GAP Algorithm 1 has been used over a grid with 2.5 km of spacing, which is the resolution of the MEPS grid (see Sec. 3). The parameters are ε 2 = 0.1, α = 0.1, 470 p mx = 200, L i = 50 km constant, D i is estimated adaptively on the grid as the distance between the grid point and the 10th closest observation location with upper and lower bounds of 3 km and 10 km, respectively. The EnSI-GAP settings are such that the analyses would stay much closer to the observations than to the forecasts, where observations are available. The analysis uncertainty will reflect locally both the forecast ensemble spread and the averaged innovation, but the weight of this last component is damped by α. The two parameters p mx and D i are used to limit the number of observations that can influence 475 the analysis at a grid point. The localization parameter L i is set to a rather large value, such that the dynamics of the forecasts ensemble are evident in the results. The observation error covariance matrix of Eq. (6) is defined with a diagonal i Γ o , that is a situation where radar-derived and in-situ observations are assumed to have the same precision, moreover we are ignoring the spatial correlation of radar-derived observation errors. An investigation of spatially correlated radar-derived observation errors is outside the scope of this study. Note that those settings are useful for the illustration of the method, while for operational 480 applications other settings may be more appropriate, such as a smaller value of L i or a more sophisticated characterization of the observation errors, for example. Figure 6 shows the hourly precipitation data for 2019-07-30 15:00 UTC over South Norway. The observational data are shown in panel a. For each grid box, the average of radar-derived precipitation and in-situ measurements within that box is Grid points that are not covered by observations are marked in gray. In panel b, the background ensemble mean derived from a 10-member ensemble forecast is shown, while six of the ten ensemble members are shown in Fig. 7. The 10-member ensemble shows realistic precipitation fields, moreover they are rather similar, at least in terms of weather situation at the Meso-β scale.
Weather forecasters can be quite confident in stating that heavy precipitation is likely to occur over western and southern Norway, while is less likely over eastern Norway. The forecast uncertainty is large enough that it is difficult to predict exactly 490 which subregion will be affected by the most intense showers. The observations confirm that showers occur along the coast of western Norway and that the most intense precipitation event is located in Sogn og Fjordane, that is the black box in Fig. 5.
Note that approximately half of the box is not covered by observations. Panel c of Fig. 6 shows the analysis, specifically the analysis mean at each grid point. In this case, the spatial analysis acts almost as a "gap filling" procedure to fill in empty spaces in between observations with the most likely precipitation values. As prescribed by our EnSI-GAP settings, the analyses over 495 observation-dense regions are not that different from the observed values.
It is interesting to have a look at the results over Sogn og Fjordane and to focus on the Meso-γ scale. The evolution in time of the hourly precipitation fields is shown in Fig. 8 for observations, background and analysis at three different hours, as clearly indicated in the figure. The two crosses mark the locations of points A and B (see Fig. 5). As stated above, point A is in a densely observed part of the domain, while point B is almost in the middle of the observation-void region and the closest 500 observations are located at a distance of approximately 10 km. As a consequence, in the EnSI-GAP algorithm, D i at point A is closer to 3 km, while at point B is closer to 10 km. The observed fields show a large variability over short distances and the difference between two adjacent points can be as large as 30 mm/h. The background is smoother than the observed field and shows scattered showers for 14:00 UTC and 15:00 UTC, then a wider precipitation cell almost centered over point B is shown at 17:00 UTC. At 14:00 UTC, the observed value at point A (from radar-derived estimates) is over 30 mm/h and it is 505 evident a sharp gradient from south-west to north-east. The gradient is so intense that the nearby points south-west of point A, only 3 km apart, shows almost no precipitation. The background indicates that a maximum of the field can occur between point A and B. The analysis matches the observations, though smoothing out their spatial variability, such that at point A the analysis value is less than 10 mm/h, that is representative of the areal average of the nearby points. A precipitation maximum of more than 30 mm/h has been reconstructed in the analysis between points A and B, that is consistent with the gradient in the 510 observations and the pattern in the background. The situation at 15:00 UTC in the observations around point A is a bit different.
The radar-estimated precipitation is again over 30 mm/h but there are several grid points in the surroundings of point A with similar values, such that the local gradient of the field is less steep and it shows a decrease of precipitation east of point A. The background also shows that it is more likely to find intense precipitation immediately to the west of point A than to the east. In general, EnSI-GAP forces the analysis to follow more closely the observations than the background and the analysis uncertainty is smaller than those of the background. As a consequence, the timing of the precipitation onset is also better represented in the analysis. At point A, the PDF of the precipitation analysis 530 between 10:00 UTC and 13:00 UTC is a Dirac's delta function and it has the value of 0 mm/h. From 14:00 UTC onward, the analysis PDF is a gamma. During this period, on average the interquartile range (i.e., the difference between the 75th and the 25th percentiles) is 20% of the analysis value, the difference between the 90th and the 10th percentiles is 38% of the analysis and the difference between the 99th and the 1st percentiles is 70% of the analysis. From 14:00 UTC to 23:00 UTC, the observed values are within the analysis envelopes shown in Fig. 9 for 50% of the hours, that is a consistent improvement compared to 535 the background. For the other 50% of the hours, the observed values lie outside the envelopes and 14:00 UTC and 19:00 UTC are the two hours when the deviations between observations and analyses are the most evident. For those two hours, the local variability of the precipitation field is extremely large, as shown in Fig. 8 for 14:00 UTC, and the observed values at point A are sort of outliers, if compared to their neighbours. The spatial analysis finds the best estimates of true values, which are areal averages, as discussed in Sec. 2.2.1 and defined in Eq. (3), with spatial supports determined by the EnSI-GAP settings. At 540 14:00 UTC and 19:00 UTC, the representativeness errors of the observations at point A are particularly large with respect to the spatial supports of the true values, such that the corresponding observations get "filtered out" by the analysis and their values are unlikely to occur according to the analysis PDF. Note that by fine-tuning the EnSI-GAP settings, the analysis PDF can be modified such that the analysis spread would become larger, which in this case would correspond to a reduction of the spatial support for the true values, and the analysis envelope would be more likely to include the observations. The trade-off between 545 e.g. accuracy and precision of the analysis at a point ultimately depends on the objective of an application. With respect to the precipitation yes/no distinction, from 14:00 UTC to 23:00 UTC, the analysis clearly shows that precipitation is occurring at the point, while the background is more uncertain. At point B, the analysis uncertainties between 10:00 UTC and 12:00 UTC is so small that the analysis PDF is a Dirac's delta function with 0 mm/h, despite there are no observations exactly located at that point. From 13:00 UTC onward, the analysis follows a gamma PDF and the spread is wider at point B that at point A. On 550 average, the interquartile range is 82% of the analysis value, the difference between the 90th and the 10th percentiles and the 99th and the 1st percentiles are both higher than 100%. However, the values depend on the weather situation, for example at the precipitation peak of 32 mm/h at 17:00 UTC, the interquartile range is 35% of the analysis value and the difference between the 90th and the 10th percentiles is 68% of the analysis. The increase in the analysis spread at point B compared to point remarkable that even for observational dense regions, such as at point A, the analysis spread remains quite large. This is due to the large variability of the observations at the Micro-scale, which by definition includes spatial scales smaller than the Meso-γ.
With respect to the EnSI-GAP theory, this large variability is interpreted as a large observation representativeness error. A very dense observational network, that is with observations that are closer than the effective resolution of the background, has two effects on EnSI-GAP: (i) it improves the accuracy of the analysis; (ii) it may increase the analysis uncertainty, such that the tails 560 are more representative of the extremes that occur in the observational dense region, possibly on scales that are not properly resolved by the analysis grid.
One of the main innovation of EnSI-GAP compared to traditional spatial analysis methods (Hofstra et al., 2008) is the specification of anisotropic background error covariances between grid points through non-stationary covariance matrices. Two visual representations of the correlations associated with those covariances are shown in Fig 10 for points A and B. The corre-565 lations are shown instead of the covariances because we are interested in the shape of the covariance patterns and correlation is a quantity which is then more correct to compare between the two points. The domain considered is Sogn og Fjordane (see Fig. 5) and the covariances are computed over the same 2.5 km grid used in Fig. 8. Then, for visualization purposes, the correlations have been downscaled over a finer resolution grid to highlight asymmetries. The rather large localization length L i = 50 km allows for the dynamics of the forecast ensemble to be visible in the analysis. The corresponding observations, In Fig. 10, the background error correlations between the points A and B and their surroundings within Sogn og Fjordane 575 are displayed. The closest two hundreds observations are shown with different symbols, depending on the rain occurrence.
The two points are only 14 km apart, nonetheless the two maps in Fig. 10 are rather different. For point A, the correlation extends westward, while it decays faster moving eastward. The area where the correlation is higher than 0.6 is confined within approximately 5 km in any direction from point A. In the corresponding analysis in Fig. 8, the shape of the area with precipitation rate higher than 30 mm/h looks like the pattern of correlation higher than 0.6 of Fig. 10. At point B, the 580 correlations tend to be more isotropic than for point A and it looks like the correlation extends more eastward than westward.
The observations 20 km north-east of point B, that are reporting no precipitation, do have correlations with point B that are comparable to those of the observations at 10 km west of it. The analysis at point B, as shown in Fig. 8, takes into account those correlations and the predicted values in the region between points A and B decrease gradually when moving from the west to the east. Instead, because of the expected better quality of those measurements, they have been reserved as independent observations for verification. This cross-validation strategy is widely used in atmospheric sciences (Wilks, 2019).

Validation over South Norway through cross-validation experiments
As for Sec. 3.1.2, the EnSI-GAP Algorithm 1 has been used. For this application, the spatial analysis predicts values at those station locations used for cross-validation. Some of the parameters are kept fixed, while others are allowed to vary. The fixed 595 parameters are: p mx = 200, L i = 50 km. Then, as in Sec. 3.1.2, D i is estimated adaptively at each location as the distance between that point and the 10th closest observation location with upper and lower bounds of 3 km and 10 km, respectively.
The parameters that are allowed to vary and that are the objective of the sensitivity analysis that follows are: ε 2 and α. There is an important difference here, in this example the radar-derived estimates are assumed to be less precise than the in-situ observations but more precise then the background. The in-situ observations are assumed to be ten times more precise than 600 the background, then ε 2 is set to 0.1 as in Sec. 3.1.2. However, the radar-derived observations are assumed to be only two times more precise than the background, or in other words they are five times less precise than the in-situ observations, and the elements of the diagonal matrix  the "yes" event for either observations or predictions is that the corresponding value must be higher than the precipitation threshold specified on the x-axis. For all predictions, it is more likely that a predicted "yes" event corresponds to an observed "yes" event for smaller thresholds than for the higher ones. The added value of the analysis over the background is evident for all configurations. The two configurations with α = 1 present similar ETS curves, though the one with ε 2 = 0.1 performs better. The same holds true when α = 0.1, though in this case the ETS is more sensitive to variations in ε 2 and the analysis 635 performance decreases faster with the increase of ε 2 .
The relative skill of the analysis probabilistic predictions over the background, used as reference, is shown by means of the Brier Skill Score (BSS) in Fig. 13. As for the ETS, also the BSS is shown as a function of the threshold used to define a "yes" event. The BSS is higher for light precipitation, then it gradually decays with the increase of the precipitation rate. The analyses with α = 1 show better performances over a wider range of amounts than those obtained with α = 0.1.

640
The reliability diagrams plotting the observed frequency against the predicted background and analysis probabilities are shown in Fig. 14 The analyses tend to cross the diagonal, which means that they are overpredicting some probabilities and underpredicting other probabilities. When α = 1, the analyses are underpredicting low probabilities and overpredicting high probabilities. The analysis configuration with α = 1 and ε 2 = 0.1 stays rather close to the diagonal, 650 overpredicting low probabilities and underpredicting high probabilities.

Discussion
For the presented implementation of EnSI-GAP, the Gaussian anamorphosis g() of Eq. (16) is based on the same gamma distribution parameters for the whole domain. This assumption might be too restrictive for very large domains, such as for all Europe for instance. In this case, different solutions may be explored such as slowly varying the gamma parameters in space, are not too large.
The EnSI-GAP implementation in Algorithm 1 requires the specification of four parameters: D, the length scale of the scale matrix; L, the localization length scale, which is governing the covariance suppression rate with the distance; α, the stabilization factor; ε 2 , used to determine the error variances. EnSI-GAP is designed such that the four parameters are location dependent. It is important to avoid abrupt variations in space of these parameters, otherwise the analysis field will show unrealistic patterns D depends on the observation spatial distribution, because it is used to ensure that the scale matrix is based on the interaction of some observations and not only just one or two. If we assume that it is reasonable to use the observational network to refine the 680 effective resolution of the background, then we can imagine that L should be larger than D. ε 2 varies between 0.1 and 1, for the reasons discussed in Sec. 2.2.2. In Sec. 3.2, we have also specified the observation errors as a function of the data source: radar data has a different quality than in-situ. It should also be possible to setup the analysis such that citizen observations have a different quality than observations measured by professional stations. We have assumed that the background ensemble is more likely to overestimate the spread than to underestimate it, for this reason we have assumed values of α between 0.1 and 1. Note that there is a further parameter that can be specified in EnSI-GAP, that is p mx , which is used to keep the analysis procedure local. However, p mx has not been considered among the four parameters optimized in Secs. 3.1.1-3.2. Instead, it has been set to 200, that is a large number for our applications. If p mx is too small, the analysis field will include discontinuities due to the sudden influence of observations that have been given significantly different weights in the analysis at a grid point with respect to the analyses at its neighbouring grid points. p mx is important for operational application over computational systems with 690 limited resources, because it can be regarded as a parameter that limits the use of computational resources within a predefined range.

Conclusions
The ensemble-based statistical interpolation with Gaussian anamorphosis (EnSI-GAP) applies inverse problem theory to the spatial analysis of hourly precipitation. Numerical model output provides the prior information, and specifically we have 695 considered ensemble forecasts, that have been combined with radar-derived estimates and in-situ observations. EnSI-GAP has been applied on datasets that are typically available within national meteorological services. In addition, opportunistic sensing networks based on citizen observations have been considered. The precipitation representation is a synthesis of all the data available. Thanks to the diffusion of open data policies, the same datasets are also nowadays available in real-time to the general public. For instance, the Norwegian Meteorological Institute provides free access to the weather forecasts and the 700 radar data used in this article via thredds.met.no, while in-situ observations, except the citizen observations, are available via frost.met.no.
EnSI-GAP assumes the precipitation fields to be locally stationary, trans-gaussian random fields. The marginal distribution of precipitation at a point is a gamma distribution. Gaussian anamorphosis is used to pre-process data in order to better comply with the requirements of linear filtering. The inverse transformation returns the gamma shape and rate. A special case is 705 considered where uncertainties are so small that the returned analysis values have delta functions as their marginal distributions.
EnSI-GAP considers each hour independently and it requires the specification of four parameters that can vary across the domain. The implementation is designed to run in parallel on a grid point by grid point basis. Despite the small number of parameters to optimize, the spatial analysis scheme is flexible enough that it can be applied also when the background ensemble is not representing the truth satisfactorily. An important case is when, in a region, all the ensemble members show 710 no precipitation, while the observations report precipitation. By adding a scale matrix to the flow-dependent background error covariance matrix, the analysis can predict precipitation even where the background is sure that it is not occurring.
The examples of applications presented allow for a better understanding of the characteristics of EnSI-GAP and they show how the statistical interpolation can be adapted to meet specific requirements. It can be used to fill in the gaps between observation-rich regions to obtain a continuous precipitation field. The analysis expected value is available everywhere, as 715 it is the background, and in observation-dense regions it can be as accurate as the observations. Thanks to the data transformation, the spread of the analysis PDF is less likely to become unrealistically large because of either large model errors or large variability of observed small-scale precipitation. Within certain limits, determined by the spatial distribution of the observa-tional network, the analysis envelope at a point can be tuned such that it is representative of the distribution of precipitation values determined by atmospheric processes occurring at smaller spatial scales than those resolved by the background. For in-720 stance, in an observation-void region, the EnSI-GAP analysis PDF at a point provides a better estimate than the background for the probability of precipitation exceeding a threshold by an observation hypothetically placed at that point. This is an important result, especially when high-impact weather is involved. Author contributions. CL developed EnSI-GAP, tested it on the case studies and prepared the manuscript with contributions from all coauthors. TN and IS configured EnSI-GAP to work with MET Norway's datasets, collected in-situ observations from opportunistic sensing networks and quality controlled them. CE prepared the radar data.   Table 2. Overview of variables and notation for local variables. All variables are specified in the transformed space. All the vectors are column vectors if not otherwise specified. If X is a matrix, Xi is its ith column (column vector) and Xi,: is its ith row (row vector).