Geostatistical interpolation by quantile kriging

The widely applied geostatistical interpolation methods of ordinary kriging (OK) or external drift kriging (EDK) interpolate the variable of interest to the unknown location, providing a linear estimator and an estimation variance as measure of uncertainty. The methods implicitly pose the assumption of Gaussianity on the observations, which is not given for many variables. The resulting “best linear and unbiased estimator” from the subsequent interpolation optimizes the mean error over many realizations for the entire spatial domain and, therefore, allows a systematic under-(over-)estimation of the variable in regions of relatively high (low) observations. In case of a variable with observed time series, the spatial marginal distributions are estimated separately for one time step after the other, and the errors from the interpolations might accumulate over time in regions of relatively extreme observations.

Therefore, we propose the interpolation method of quantile kriging (QK) with a two-step procedure prior to interpolation: we firstly estimate distributions of the variable over time at the observation locations and then estimate the marginal distributions over space for every given time step. For this purpose, a distribution function is selected and fitted to the observed time series at every observation location, thus converting the variable into quantiles and defining parameters. At a given time step, the quantiles from all observation locations are then transformed into a Gaussian-distributed variable by a 2-fold quantile–quantile transformation with the beta- and normal-distribution function. The spatio-temporal description of the proposed method accommodates skewed marginal distributions and resolves the spatial non-stationarity of the original variable. The Gaussian-distributed variable and the distribution parameters are now interpolated by OK and EDK. At the unknown location, the resulting outcomes are reconverted back into the estimator and the estimation variance of the original variable. As a summary, QK newly incorporates information from the temporal axis for its spatial marginal distribution and subsequent interpolation and, therefore, could be interpreted as a space–time version of probability kriging.

In this study, QK is applied for the variable of observed monthly precipitation from raingauges in South Africa. The estimators and estimation variances from the interpolation are compared to the respective outcomes from OK and EDK. The cross-validations show that QK improves the estimator and the estimation variance for most of the selected objective functions. QK further enables the reduction of the temporal bias at locations of extreme observations. The performance of QK, however, declines when many zero-value observations are present in the input data. It is further revealed that QK relates the magnitude of its estimator with the magnitude of the respective estimation variance as opposed to the traditional methods of OK and EDK, whose estimation variances do only depend on the spatial configuration of the observation locations and the model settings.

1 Introduction

Many environmental variables (e.g., precipitation, ore grades) are only measured at some distinct observation locations, but possess a highly variable and unknown spatial distribution (Armstrong, 1998). The regionalization of the variable, i.e., the interpolation of the variable to the unknown locations, attempts the full description of its spatial distribution as a prerequisite for practical objectives (e.g., hydrological modeling, efficient exploitation of resources). However, a deterministic description of the spatial distribution is severely hampered for many variables since they incorporate a complex genesis which is neither fully known nor understood.

Therefore, the assessment of the distribution by geostatistical models arose, whose theoretical fundamentals were firstly laid out by Matheron (1962). The theory regards the (observed) values of a variable z as one realization z(x) of a random variable (RV) Z at the specific location x $(=\mathit{\{}{x}_{\mathrm{1}},x_{\mathrm{2}}\mathit{\}}$ for ℝ²). Since the variable is often only partially known at a few distinct measurement locations x_i, ergodicity is assumed and the intrinsic hypothesis for ordinary kriging (OK, Matheron, 1965) is given as

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }E\left[Z\right(\mathit{x}\left)\right]=m\left(\mathit{x}\right)=\mathrm{const}.,\text{(2)}& {\displaystyle }& {\displaystyle }\mathrm{VAR}\left[Z\right(\mathit{x}+\mathit{h})-Z(\mathit{x}\left)\right]=\mathrm{2}\mathit{\gamma }\left(\mathit{h}\right),\end{array}$$

in which γ(h) is the semi-variogram. The increment $Z(\mathit{x}+\mathit{h})-Z\left(\mathit{x}\right)$ is assumed as a stationary random function and its variance only depends on the translation vector h. The set of outcomes by the interpolation to the unknown location x is described by the linear estimator Z^∗(x) as a measure of their centrality and by the estimation (or kriging) variance ${\mathit{\sigma }}_{\mathrm{K}}^{\mathrm{2}}\left(\mathit{x}\right)$ as a measure of the associated uncertainty.

The stated hypothesis entails three implications: the first condition of the intrinsic hypothesis (Eq. 1) demands the variable to be spatially stationary and yields an unbiased estimation error in the entire domain (Chilès and Delfiner, 1999). Therefore, a systematic under-(over-)estimation is induced for regions of high (low) observations. Secondly, the marginal distribution of the observed data should ideally be Gaussian in order to be adequately described. Unfortunately, the distribution often departs from the ideal case (Journal and Alabert, 1989), necessitating an a priori transformation of the marginal distribution. And last, the second condition (Eq. 2) implies that the magnitude of the kriging variance ${\mathit{\sigma }}_{\mathrm{K}}^{\mathrm{2}}\left(\mathit{x}\right)$ only depends on the spatial configuration of the observation locations, the a priori variance of all observations and on the selected variogram model, but not on the magnitude of the linear estimator Z^∗(x) itself (Goovaerts, 2000).

The theoretical extension of external drift kriging (EDK, Ahmed and deMarsily, 1987) addresses the first implication. EDK can be attributed to the non-stationary geostatistical interpolation methods and has been frequently applied in various disciplines of practice and science (e.g., Bourennane et al., 2000, van de Kassteele et al., 2009, Motaghian and Mohammadi, 2011). It incorporates additional information from external variables (or drifts) Y_j(x) for the estimation of the variable at the unknown location. The mean m(x) is non-stationary but linearly depends on the external variables. Thus, the first condition (Eq. 1) is reformulated to

$$\begin{array}{}\text{(3)}& E\left[Z\left(\mathit{x}\right)\right]=a+\sum _{i=\mathrm{1}}^k{b}_i\cdot Y_i\left(\right(\mathit{x}\left)\right),\end{array}$$

where k is the number of the incorporated external drifts Y_i(x), while a, b_i are the unknown constants. The drifts are required to be known prior to interpolation at all relevant locations. Ideally, EDK still requires a marginal Gaussian distribution.

The non-parametric methods of indicator kriging and probability kriging (Journel, 1983) are further derivatives and allow an a priori transformation of the skewed marginal distribution. Indicator kriging transforms the variable z into a binary variable by defining a threshold value z_th and restates the first condition (Eq. 1) of the random function $I(\mathit{x};z_{\mathrm{th}})={\mathrm{1}}_{Z\left(\mathit{x}\right)<z_{\mathrm{th}}}$ to

$$\begin{array}{}\text{(4)}& E\left[I(\mathit{x};z_{\mathrm{th}})\right]=m=\mathrm{const}.\end{array}$$

This nonlinear transformation is relatively robust and limits the effect of high values on the description of the variable at the unknown location. However, a loss of information comes along by the transformation into a binary variable. Therefore, probability kriging defines multiple thresholds and uses the order relation of the observed variable (Carr and Mao, 1993), being implemented by using co-kriging in the derived multi-variate context. Both non-parametric methods have been subject to research, especially for detection limit problems of groundwater contamination (e.g., Goovaerts et al., 2005; Lee et al., 2007; Adhikary et al., 2011).

In summary, geostatistical methods have been derived in the past in order to address the stated shortfalls of the intrinsic hypothesis. However, all present methods only regard the observations from the one respective time step for the estimation of their marginal spatial distribution, but do not incorporate observations from other time steps. The inclusion of a temporal behavior into the geostatistical models is mostly irrelevant for the original geological variables. However, the temporal variability of a variable becomes more prominent for other sciences, e.g., hydrology, where observations from raingauges over several time steps are implemented into the geostatistical models in order to generate spatial precipitation estimates. These estimates subsequently serve as input to the hydrological modeling (e.g., Syed et al., 2003; Basistha et al., 2008; Cole and Moore, 2008) over multiple time steps. Associated errors in the precipitation estimates may ultimately lead to greater errors in the subsequent discharge modeling (Kobold and Sušelj, 2005). These errors strongly depend on the spatial and temporal distribution of the input precipitation (Gabellani et al., 2007; Moulin et al., 2009) and may limit the accuracy of rainfall–runoff simulations. There are geostatistical space–time models in order to incorporate the temporal variability of the variable, but they are primarily aiming on the extrapolation of the variable in time (Snepvangers et al., 2003). Therefore, they require a strong dependence of the variable over time, suited, e.g., for groundwater modeling where temporal changes occur relative slowly. This temporal dependence might be absent for other variables, e.g., monthly precipitation.

In the following section, we introduce quantile kriging as a spatio-temporal description of the variable Z, addressing the three shortfalls: non-stationary variables, skewed marginal distributions over space and the independence of the error distribution from the magnitude of the observation.

2 Materials and methods

The theory of quantile kriging (QK) is outlined along with the major theoretical implications, followed by a general discussion of the underlying geostatistical model and a case study for the variable of monthly precipitation is presented.

Figure 1 Histogram from the times series of observed monthly precipitation for two random stations: “Laingsnek”, x₅₇ with $\stackrel{\mathrm{‾}}{z}=\mathrm{81.3}$ mm (a) and “Tambotieboom”, x₂₉ with $\stackrel{\mathrm{‾}}{z}=\mathrm{38.3}$ mm (b).

A preliminary analysis of the selected variable exemplary reveals (Fig. 1) the non-Gaussianity within the data and that the first assumption of the intrinsic hypothesis (see Eq. 1) is not fulfilled since, e.g., $E_t\left[Z\right({\mathit{x}}_{\mathrm{57}},t\left)\right]\ne E_t\left[Z\right({\mathit{x}}_{\mathrm{29}},t\left)\right]$.

2.1 Theory of quantile kriging

QK presumes the existence of observations of the variable z over consecutive time steps t $(=\mathrm{1},\mathrm{2},\mathrm{\dots },J)$ at every observation location x_i $=(x{\mathrm{1}}_i,x{\mathrm{2}}_i)$ (for the two-dimensional space ℝ²), providing an observed time series z(x_i,t) at every observation location i $(=\mathrm{1},\mathrm{2},\mathrm{\dots },n_i)$. QK proceeds first with a two-step procedure prior to interpolation (see Sect. 2.1.1) and second with the interpolation itself (see Sect. 2.1.2).

2.1.1 Estimation of the temporal and the spatial marginal distribution

At first, the distribution over time is estimated at every observation location location x_i: an appropriate theoretical cumulative distribution function (cdf) F is selected and fitted to the corresponding time series of observations z(x_i,t), yielding n_i specific distributions F(z(x_i,t)). The distributions are defined by their corresponding parameter sets Θ(x_i) (= ϑ_k(x_i) with $k=\mathrm{1},\mathrm{2},\mathrm{\dots },K$) of the K-parametric distribution function F. The quantiles w(x_i,t) (= $F\left(z\right({\mathit{x}}_i,t);\mathrm{\Theta }({\mathit{x}}_i)$) are calculated from the observed values of the variable z(x_i,t) and the defining parameter set Θ(x_i). The quantiles w(x_i,t) possess a uniform distribution over time on the interval [0,1] for a given observation location x_i. However, their empirical distribution in space is not uniform on [0,1]. In order to profit from the optimality of kriging, it requires a transformation into a Gaussian distribution as a prerequisite of the subsequent interpolation.

The marginal spatial distribution corresponding to a time step t is, therefore, estimated by a 2-fold quantile–quantile conversion as the second step: the two-parametric beta distribution is fitted to the quantiles w(x_i,t) of a given time step t, whose cdf $G(w;\mathit{\alpha },\mathit{\beta })$ is defined as

$$\begin{array}{ll}\text{(5)}& {\displaystyle }G(w;\mathit{\alpha },\mathit{\beta })& {\displaystyle }=\sum _{n=\mathit{\alpha }}^{\mathit{\alpha }+\mathit{\beta }-\mathrm{1}}{\displaystyle }& {\displaystyle }\left[{\displaystyle \frac{(\mathit{\alpha }+\mathit{\beta }-\mathrm{1})\mathrm{!}}{n\mathrm{!}\cdot (\mathit{\alpha }+\mathit{\beta }-\mathrm{1}-n)\mathrm{!}}}\cdot w^n\cdot (\mathrm{1}-w{)}^{(\mathit{\alpha }+\mathit{\beta }-\mathrm{1}-n)}\right]\end{array}$$

on the interval [0,1] by the two parameters α>0 and β>0. The quantiles $G\left(w\right({\mathit{x}}_i,t);\mathit{\alpha }(t),\mathit{\beta }(t\left)\right)$ from Eq. (5) are finally transformed by a normal score transformation into the standard Gaussian variable u(x_i,t) with N_u[0|1], which ultimately serves as spatial marginal distribution to the subsequent geostatistical interpolation.

The transformation via the quantiles into the variable u accounts for spatially non-stationary distributions of the original variable z with $E\left[Z\right(\mathit{x},t\left)\right]\ne m$ and exchanges the two conditions of Eqs. (1) and (2) to

$$\begin{array}{}\text{(6)}& {\displaystyle }& {\displaystyle }E\left[F_x\right(Z(\mathit{x},t)\left)\right]=m=\mathrm{const}.,\text{(7)}& {\displaystyle }& {\displaystyle }\mathrm{VAR}\left[F_{x+h}\right(Z(\mathit{x}+\mathit{h},t))-F_x(Z(\mathit{x},t)\left)\right]=\mathrm{2}\mathit{\gamma }\left(\mathit{h}\right),\end{array}$$

resolving the problem of spatial non-stationarity. QK can accommodate skewed marginal distributions of the original variable Z, which is similar to probability kriging, but it newly incorporates the temporal behavior of Z into its estimation of the spatial marginal distribution.

2.1.2 Interpolation to the unknown location

The outlined conversion of the variable z(x_i,t) into the variable u(x_i,t) and its corresponding parameter set Θ(x_i) entails separate interpolations to the unknown location x.

The inherent assumption of second-order stationarity implies the existence of a constant spatial mean for the variable u within the domain for every time step t. The transformed quantiles are implicitly assumed to be more homogeneously distributed over space than the original variable z. The variable u is subject to a stationary geostatistical interpolation method (e.g., OK), providing a linear estimator $U^{\ast }(\mathit{x},t)$ and the estimation variance ${\mathit{\sigma }}_{K,U}^{\mathrm{2}}(\mathit{x},t)$. They jointly describe the Gaussian distribution of the random variable U(x,t) with $N\left[U^{\ast }\right|{\mathit{\sigma }}_{K,U}^{\mathrm{2}}]$.

The defining parameters ϑ_k (for $k=\mathrm{1},\mathrm{2},\mathrm{\dots },K$) of the K-parametric distribution function F are independent from the time step t and they are separately interpolated to the unknown location x. The separate interpolation, however, requires the independence of the parameters ϑ_k from each other. Therefore, a principal component analysis examines the “observed” parameters ϑ_k(x_i) in the Cartesian coordinate system (${\mathit{\vartheta }}_{\mathrm{1}}\left|{\mathit{\vartheta }}_{\mathrm{2}}\right|\mathrm{\dots }|{\mathit{\vartheta }}_{\mathrm{K}}$) and determines the corresponding rotation angle α and translation vector k. The coordinate system is then rotated and translated accordingly prior to interpolation in order to ensure independence. A possible spatial non-stationarity of the parameters can be accounted for by the choice of an appropriate non-stationary interpolation method (e.g., EDK). The interpolation of the independent parameters yields their estimators at the unknown location x, which are rotated and translated back to the original coordinate system. Thus, the estimators ${\mathit{\vartheta }}_k^{\ast }\left(\mathit{x}\right)$ are defining the distribution function at the unknown location x.

At last, the resulting Gaussian distribution of the random variable U(x,t) is reconverted into a distribution of the original variable Z(x,t) by the outlined steps of conversion (Sect. 2.1.1), but in reverse order: first the distribution of the quantiles $G\left(W\right(\mathit{x},t);\mathit{\alpha },\mathit{\beta })$ of the beta distribution is calculated by the inverse of the normal score transformation. The distribution of the quantiles W(x,t) are calculated next using the inverse of Eq. (5) and last, the distribution of the original variable Z(x,t) is estimated by using the inverse of the selected cdf, being defined by the estimators of its parameters ${\mathit{\vartheta }}_k^{\ast }\left(\mathit{x}\right)$. The reconversion of the distribution of U(x,t) to the distribution of Z(x,t) can be implemented by the simple numerical Rosenblueth point estimation method (Rosenblueth, 1975). The resulting distribution of the original variable Z is then described by the expectation value $Z^{\ast }(\mathit{x},t)$ and the variance ${\mathit{\sigma }}_{\mathrm{K}}^{\mathrm{2}}(\mathit{x},t)$. Note that the resulting asymmetrical distribution of Z(x,t) is non-Gaussian due to the conversion with the nonlinear but monotonic theoretical cdf F.

Figure 2 Flowchart for the basic methodology of quantile kriging.

The basic methodology of the proposed QK is illustrated in Fig. 2.

2.2 Discussion of the geostatistical model

Since the proposed method of QK is applied for the variable of monthly precipitation (see Sect. 2.3), the discussion of the underlying process model is based on the following properties of precipitation fields:

• The monthly (and even daily) precipitation amounts z(x,t) for a given time step t often show a skewed distribution and cannot be considered as stationary over space. The differences in expected precipitation amounts become especially obvious for long time accumulations.

• The meteorological processes, which are generating precipitation, are usually of large spatial extent: if one location receives heavy precipitation, it is likely that other locations also receive heavy precipitation.

• Correlations between time series of precipitation indicate a strong spatial dependence, while the spatial dependence of precipitation at one given time step (e.g., day, month) usually show a much weaker spatial dependence.

A possible process model reflecting the above properties can be described as follows.

Let U(x,t) be independent (for each different time step t) normal stationary spatial fields with E[U]=0 and D²(U)=1 for each time step t. Now, the process M(t) is introduced in order to reflect large-scale meteorological processes. High M(t) values correspond to heavy precipitation covering the area, while low values correspond to dry conditions, as it is reflected by seasonal variations of precipitation amounts. The introduced M modifies the spatial process to

$$\begin{array}{}\text{(8)}& G(\mathit{x},t)=U(\mathit{x},t)+M\left(t\right),\end{array}$$

where M(t) is a process (only in time) with a mean of zero. We may assume that the distribution of M(t) is normal and, therefore, G(x,t) would be normally distributed with N(0,d) (with $d=\sqrt{\mathrm{1}+{\mathit{\sigma }}_M^{\mathrm{2}}}$) at every location x.

For each individual time step t, the distribution of G(x,t) is N(M(t),1) and the resulting spatial field W is the temporal non-exceedance probability at location x being confined to $\mathrm{0}\le W(\mathit{x},t)\le \mathrm{1}$ and formally described as

$$\begin{array}{}\text{(9)}& W(\mathit{x},t)={\mathrm{\Phi }}_{\mathrm{0},\mathrm{d}}\left(G\right(\mathit{x},t\left)\right),\end{array}$$

where Φ_0,d is the distribution function of N(0,d). The precipitation is then generated as

$$\begin{array}{}\text{(10)}& Z(\mathit{x},t)=F_x^{-\mathrm{1}}\left(W\right(\mathit{x},t\left)\right),\end{array}$$

where F_x is the distribution function of precipitation at the location x. The distribution functions F_x may vary between different locations x due to topography and other influencing factors, and they could be subject to interpolation (e.g., Mosthaf and Bárdossy, 2017).

We use W(x,t) for each time step t and assume that it follows a beta distribution. In fact, its distribution depends on M(t). If M(t)=0 for all time steps t, then monthly precipitation can be fully characterized by independent realizations over space. In this case, the distribution of W is uniform for each t.

However, this is not the case with observed data because wet and dry conditions occur simultaneously over the entire domain. This is controlled by M(t), which can be taken, e.g., as an independent random variable or to follow an ARMA process. If M(t)≠0 then the distribution of W(x,t) is not uniform for this specific time step t. The exact form of the corresponding distribution would be something like

$$\begin{array}{}\text{(11)}& G_t\left(v\right)={\mathrm{\Phi }}_{\mathrm{0},\mathrm{1}}\left({\mathrm{\Phi }}_{M\left(t\right),\mathrm{1}}^{-\mathrm{1}}\left(v\right)\right).\end{array}$$

However, the use of Eq. (11) would require the estimation of M(t) for each time step t. We decided to use a simple beta distribution instead. The reason for assuming a beta distribution is due to their flexibility and their ability to describe distributions well within the interval [0,1].

The introduction of M(t) is reasonable as it explains the difference between the correlation between stations and the spatial correlation calculated using a variogram type approach for a given time step. The later correlations are usually lower (smaller ranges), which are increased by the common large-scale weather described by M(t). Note that the introduction of M(t) leads to a correlation of the precipitation time series even if the individual snapshots of U(x,t) are independent in space.

We estimate and subsequently interpolate F_x within the proposed methodology by the preceding conversion of the variable Z(x_i,t). In addition, we calculate W(x_i,t) for the observation locations x_i and interpolate it to the unknown location x in order to come back to Z(x,t). In here, spatial variograms are calculated for W for each time step t, assuming W to be spatially stationary. Non-Gaussian and non-stationary distributions only occur for the precipitation amounts (i.e., the variable Z).

Non-Gaussianity should be considered due to the usually skewed distribution of precipitation amounts and it only applies to the marginal distribution at a given time step t. The suggested model should enable a simulation of the precipitation amounts. The spatial dependencies are considered to correspond to a multi-Gaussian copula, being a type of transformation frequently used (e.g., for lognormal kriging).

The distributions F_x, fitted to the individual locations, are supposed to have a spatial dependence. They are further assumed to follow the same distribution (e.g., Γ or Weibull distribution) and are subsequently interpolated. In here, we assume that the large-scale meteorological processes, generating precipitation, are better reflected by the distributions than by a single monthly (or daily) realization. Therefore, the use of external covariates, e.g., elevation, is deemed more appropriate for their interpolation. The usage of these distributions transforms the process into a stationary one, which is then interpolated using the beta distribution of the non-exceedance probabilities.

2.3 Application of quantile kriging

The proposed method of QK is applied for the variable of monthly precipitation in South Africa and the outcomes are compared to those from OK and EDK.

2.3.1 Study area and data

The rectangular study area ($\mathrm{3.5}{}^{\circ }\times \mathrm{3.5}{}^{\circ }$, Fig. 3) covers approx. 132 000 km² and is located within the Republic of South Africa. The second release of the digital elevation model from the Shuttle Radar Topography Mission (USGS, 2003) serves as elevation input. The original resolution was upscaled from 3 arcsec (approx. 92 m) to 2 arcmin (approx. 3700 m) by spatial averaging, resulting in a mean of 1442 m and ranging from 669 to 2197 m (a.m.s.l.). The upscaled elevation ultimately serves as external drift for EDK of the parameters within QK and for the reference EDK with the original variable.

Figure 3 Study area, elevation and location of raingauges.

The observations of monthly precipitation were retrieved from raingauges of four different sources: the Department of Water Affairs (DWA, 2008), the Global Historical Climatology Network (Vose et al., 1992), the Climate Research Unit (Mitchell and Jones, 2005) and the internal database of the University of KwaZulu-Natal (Lynch, 2004). Accumulation of daily recordings yield monthly values for the 264 (=J) consecutive months from January 1986 to December 2007. A total of 226 (=n_i) raingauges (Fig. 3) provided 32 226 (=N) monthly precipitation values, which ultimately serve as input data.

The observed average monthly precipitation over the 12 calendar months c is illustrated in Fig. 4 along with the percentage of zero-value observations over all observations of the specific calender month c (hereafter referred to as the dry ratio), revealing a seasonal variation. High precipitation is typically encountered in the calendar months from October to March, being characterized by a low dry ratio <3 %. The study area receives relatively low precipitation amounts during the calender months from April (dry ratio = 11 %) to September (dry ratio = 25 %).

2.3.2 Adaptation to monthly precipitation

At first, we subdivided the observations of monthly precipitation into the corresponding calendar month c $(=\mathrm{1},\mathrm{2},\mathrm{\dots },\mathrm{12})$ prior to the fitting of the selected distribution function due to two reasons: the seasonal variation in monthly precipitation (Fig. 4) and to ensure independence of the individual sample members as a theoretical requirement for the fitting method. We used the maximum likelihood estimation method for fitting the selected distribution function to the respective measurements values z(x_i,t_c) of every calender month c and every measurement location x_i, resulting in a total of 2712 ($=\mathrm{12}\times \mathrm{226}$) fittings. In this context, the two-parametric Γ and Weibull distribution were selected, whose cdf F(z;Θ) are defined as

$$\begin{array}{}\text{(12)}& {\displaystyle }& {\displaystyle }\mathrm{\Gamma }-\mathrm{distribution}:F(z;\mathrm{\Theta })={\displaystyle \frac{\mathit{\gamma }(\mathit{\mu },\mathit{\lambda }\cdot z)}{\mathrm{\Gamma }\left(\mathit{\mu }\right)}},\text{(13)}& {\displaystyle }& {\displaystyle }\mathrm{Weibull}-\mathrm{distribution}:F(z;\mathrm{\Theta })=\mathrm{1}-\exp \left[-{\left({\displaystyle \frac{z}{\mathit{\lambda }}}\right)}^k\right],\end{array}$$

where Γ(μ) is the gamma function and $\mathit{\gamma }(\mathit{\mu },\mathit{\lambda }\cdot z)$ is the lower incomplete gamma function. The parameter set Θ_c(x_i) is composed for the Γ distribution out of ${\mathit{\mu }}_c\left({\mathit{x}}_i\right)(={\mathit{\vartheta }}_{\mathrm{1},c}({\mathit{x}}_i\left)\right)$ and ${\mathit{\lambda }}_c\left({\mathit{x}}_i\right)(={\mathit{\vartheta }}_{\mathrm{2},c}({\mathit{x}}_i\left)\right)$ and for the Weibull distribution out of $k_c\left({\mathit{x}}_i\right)(={\mathit{\vartheta }}_{\mathrm{1},c}({\mathit{x}}_i\left)\right)$ and ${\mathit{\lambda }}_c\left({\mathit{x}}_i\right)(={\mathit{\vartheta }}_{\mathrm{2},c}({\mathit{x}}_i\left)\right)$. All parameters are restrained to values greater than zero and both cdfs are defined for $z({\mathit{x}}_i,t)\ge \mathrm{0}$.

Figure 4 Average monthly precipitation (in mm) and dry ratio (in %) from 226 raingauges. Note that the dry ratio (dashed brown line) is indicated on the right axis.

Thus, the original observations of monthly precipitation z(x_i,t) are converted by Eq. (12) or Eq. (13) into the corresponding quantiles w(x_i,t) (= $F\left(z\right({\mathit{x}}_i,t);{\mathrm{\Theta }}_c)$) and their defining parameter set Θ_c(x_i). As outlined in Sect. 2.1, the quantiles w(x_i,t) were further converted into the standard Gaussian variable u(x_i,t), ultimately subject to the subsequent OK as our chosen geostatistical interpolation method. Note that the stationary assumption of more homogeneously distributed quantiles in space appear more plausible in the case of monthly precipitation. In total, the variable u(x_i,t) was interpolated 264 times by OK to the unknown location x. The corresponding variograms are calculated using Kendall's tau for a robust interpolation (Lebrenz and Bárdossy, 2017).

We further selected EDK as non-stationary interpolation method for the defining parameters ϑ_1,c and ϑ_2,c. Elevation data (Sect. 2.3.1) are taken as external drift since the distributions of monthly precipitation are assumed to depend on the altitude of the terrain. However, it was revealed that the direct use of the parameters may lead to negative or zero estimators at locations of extreme external drifts. Therefore, the sample mean ${\stackrel{\mathrm{‾}}{z}}_c\left({\mathit{x}}_i\right)$ and the sample variance ${\mathit{\sigma }}_c^{\mathrm{2}}\left({\mathit{x}}_i\right)$ are estimated instead, using the two statistical moments of the selected distribution functions, defined as

$$\begin{array}{ll}\text{(14)}& {\displaystyle }& {\displaystyle }\mathrm{\Gamma }-\mathrm{distribution}:E\left[Z\right]={\displaystyle \frac{\mathit{\mu }}{\mathit{\lambda }}};\mathrm{VAR}\left[Z\right]={\displaystyle \frac{\mathit{\mu }}{{\mathit{\lambda }}^{\mathrm{2}}}};{\displaystyle }& {\displaystyle }\mathrm{Weibull}-\mathrm{distribution}:E\left[Z\right]=\mathit{\lambda }\cdot \mathrm{\Gamma }(\mathrm{1}+\mathrm{1}/k);\\ \text{(15)}& {\displaystyle }& {\displaystyle }\mathrm{VAR}\left[Z\right]={\mathit{\lambda }}^{\mathrm{2}}\cdot \mathrm{\Gamma }(\mathrm{1}+\mathrm{2}/k)-E[Z{]}^{\mathrm{2}}.\end{array}$$

The dependence of the two derived parameters $\stackrel{\mathrm{‾}}{z}$ and σ² on each other appears obvious in the case of monthly precipitation: a high mean is likely to be associated with a high variance and vice versa. Their dependence is exemplary illustrated for the calender month “May” in Fig. 5.

Figure 5 Scatterplot of sample mean ${\stackrel{\mathrm{‾}}{z}}_{\mathrm{5}}$ and sample variance ${\mathit{\sigma }}_{\mathrm{5}}^{\mathrm{2}}$ of calendar month “May” along with the principle components s₅ and r₅.

The principal component analysis allows for the transformation into the new Cartesian coordinate system with the new coordinates r_c(x_i) and s_c(x_i). They are now independent and subject to a separate interpolation by EDK. A total of 24 ($=\mathrm{12}\times \mathrm{2}$) interpolations by EDK to the unknown location x is performed for each selected type of distribution.

3 Results and discussion

The proposed interpolation method of QK, using either a Γ distribution (QK-Γ) or a Weibull distribution (QK-Wei), is implemented and compared to the traditional geostatistical interpolation methods of OK and EDK. The respective performances are evaluated by cross-validation for the resulting estimators Z^∗ and the associated kriging variances ${\mathit{\sigma }}_{\mathrm{K}}^{\mathrm{2}}$. In here, cross-validation eliminates all values z(x_i,t) in turns from the input data, and subsequently calculates the estimator Z^∗(x_i) and the associated kriging variance ${\mathit{\sigma }}_{\mathrm{K}}^{\mathrm{2}}\left({\mathit{x}}_i\right)$ from the remaining data. Only the 32 226 data points of the actually recorded values were considered for the cross-validation and the resulting outcomes are compared to the actual observations.

3.1 Implementation of quantile kriging

The outcomes from the interpolation by OK, EDK and QK-Γ are exemplary displayed and examined for a month with low precipitation and a high dry ratio (August 1993) and a month with high precipitation and a low dry ratio (January 1996). The respective spatial patterns of the estimator Z^∗(x) and the associated standard deviation σ_K(x) are illustrated in Figs. 6 and 7.

Figure 6 Spatial patterns of the estimator Z^∗(x) and the standard deviation σ_K(x) from OK, EDK and QK-Γ for August 1993. Note that crosses indicate positions of raingauges.

The estimator Z^∗ displays similar spatial patterns and value ranges for all the interpolation methods. However, the local contours of the isohyets are more rugged for QK-Γ (Figs. 6e and 7e) than for OK (Figs. 6a and 7a), but smoother than for EDK (Figs. 6c and 7c).

QK utilizes elevation for the interpolation of the two distribution parameters ϑ_1,c and ϑ_2,c. The two parameters incorporate information from all time steps t_c of the specific calendar month c and, thus, transfer information over time. They are further combined with the ordinary kriged quantiles W(x,t), leading to more smooth contours of the isohyets than EDK (compare Fig. 3). We regard the resulting spatial patterns of QK as more plausible, assuming that the accumulated monthly precipitation is hardly affected by local features in elevation.

Figure 7 Spatial patterns of the estimator Z^∗(x) and the standard deviation σ_K(x) from OK, EDK and QK-Γ for January 1996. Note that crosses indicate positions of raingauges.

The standard deviations σ_K of the associated estimation error show notable deviations in spatial pattern for the implemented interpolation methods. The range of error is notably higher for QK (Figs. 6f and 7f) and its spatial patterns deviates from the typical, bull-eye-shaped patterns of OK (Figs. 6b and 7b) or EDK (Figs. 6d and 7d).

The estimation error from OK and EDK depends on the spatial configuration of the observation locations x_i, their global variance and the selected variogram model. This typical spatial pattern of the error distribution from the ordinary kriged quantiles W(x) is converted within QK by the monotonic cdf (Eq. 12). The resulting σ_K(x) of the original variable Z(x) is, therefore, increased by relatively flat slopes of the cdfs, which are encountered for relatively high values of W(x). A relationship between the magnitude of the estimator Z^∗ and the magnitude of the associated standard deviation σ_K is suggested by Figs. 6 and 7.

3.1.1 Relationship between estimator and standard deviation

A relationship between the magnitude of the estimator Z^∗ and the associated standard deviation σ_K would be possibly an improvement to geostatistical interpolation and is, therefore, examined next by cross-validation. The Spearman rank correlation coefficient ρ_S is chosen for its description, being defined as

$$\begin{array}{ll}\text{(16)}& {\displaystyle }& {\displaystyle }{\mathit{\rho }}_{\mathrm{S}}={\displaystyle }& {\displaystyle \frac{\sum _{i=\mathrm{1}}^n\left(\mathrm{rg}\right(Z^{\ast }({\mathit{x}}_i,t))-{\stackrel{\mathrm{‾}}{\mathrm{rg}}}_{Z^{\ast }})\times \left(\mathrm{rg}\right({\mathit{\sigma }}_{\mathrm{K}}({\mathit{x}}_i,t))-{\stackrel{\mathrm{‾}}{\mathrm{rg}}}_{\mathit{\sigma }})}{\sqrt{\sum _{i=\mathrm{1}}^n\left(\mathrm{rg}\right(Z^{\ast }({\mathit{x}}_i,t))-{\stackrel{\mathrm{‾}}{\mathrm{rg}}}_{Z^{\ast }}{)}^{\mathrm{2}}\times \sum _{i=\mathrm{1}}^n(\mathrm{rg}\left({\mathit{\sigma }}_{\mathrm{K}}\right({\mathit{x}}_i,t\left)\right)-{\stackrel{\mathrm{‾}}{\mathrm{rg}}}_{\mathit{\sigma }}{)}^{\mathrm{2}}}}},\end{array}$$

where $\mathrm{rg}\left(Z^{\ast }\right({\mathit{x}}_i,t\left)\right)$ and rg(σ_K(x_i,t)) are the ranks of the estimator Z^∗ and the associated standard deviation σ_K within a set of data, while ${\stackrel{\mathrm{‾}}{\mathrm{rg}}}_{Z^{\ast }}$ and ${\stackrel{\mathrm{‾}}{\mathrm{rg}}}_{\mathit{\sigma }}$ are the respective average ranks. The non-parametric Spearman rank correlation ρ_S describes the monotonic relation between the estimator Z^∗ and estimation standard deviation σ_K, ranging from −1 (negative) to +1 (positive) with 0 indicating its absence. A set of data consists of all n values of the corresponding calendar month c. The evolution of the Spearman rank correlation coefficient ρ_S over all 12 calendar months is displayed in Fig. 8.

Figure 8 Evolution of the Spearman rank correlation coefficient ρ_S between the estimator Z^∗ and the standard deviation σ_K.

The rank correlation varies over the calendar months for all implemented interpolation methods and reach their seasonal maximums in June or July (Fig. 8), being characterized by a high dry ratio and low precipitation.

An improvement in the relationship between the estimator Z^∗ and the associated standard deviation σ_K can be observed for QK-Γ and QK-Wei, exhibiting superior rank correlation coefficients for all calender months with the exception of QK-Wei in December (Fig. 8). QK-Γ deploys the strongest relation during the wetter months from October to March, while QK-Wei is superior from May to September. The resulting spread of the error distribution is increased by decreasing slopes of the theoretical cdfs (Eqs. 12 and 13) and vice versa. The slope is effectively the probability density function (pdf). Both selected theoretical distributions imply a monotonic decrease in their respective pdfs for small parameters, being typically encountered during the dry season, and evoke a higher spread of the error distribution for higher monthly precipitation. Thus, the almost perfect rank correlation ρ_S(c) of QK during the months of low precipitation can be explained. The rank correlation between the estimator and the standard deviation is weakened for the months with higher precipitation due to the departure from the strict monotonic decrease in the pdfs, which is induced by increasing distribution parameters.

The inferior correlation coefficients of OK and EDK are nearly congruent due to their inherent geostatistical definition: although the kriging weights are altered by the drift, they influence the linear estimator Z^∗ and the standard deviation σ_K by the same extent. Therefore, the non-parametric Spearman descriptor hardly differentiates between OK and EDK.

3.2 Cross-validation of the estimator

The estimator $Z^{\ast }({\mathit{x}}_i,t_j)$ from the cross-validation is evaluated by six objective functions: the Pearson correlation coefficient ρ, the Nash–Sutcliffe efficiency coefficient NSE (Nash and Sutcliffe, 1970), the overall bias B1 and the root mean square error (RMSE) are complemented by the temporal bias B2 and the spatial bias B3 (Bárdossy and Pegram, 2012), which are defined as

$$\begin{array}{}\text{(17)}& {\displaystyle }& {\displaystyle }\mathit{\rho }={\displaystyle \frac{\sum _{i=\mathrm{1}}^{n_i}\sum _{j=\mathrm{1}}^J\left(Z^{\ast }\right({\mathit{x}}_i,t_j)-{\stackrel{\mathrm{‾}}{Z}}^{\ast })\times \left(z\right({\mathit{x}}_i,t_j)-\stackrel{\mathrm{‾}}{z})}{\sqrt{\sum _{i=\mathrm{1}}^{n_i}\sum _{j=\mathrm{1}}^J{\left[Z^{\ast }({\mathit{x}}_i,t_j)-{\stackrel{\mathrm{‾}}{Z}}^{\ast }\right]}^{\mathrm{2}}\times \sum _{i=\mathrm{1}}^{n_i}\sum _{j=\mathrm{1}}^J{\left[z({\mathit{x}}_i,t_j)-\stackrel{\mathrm{‾}}{z}\right]}^{\mathrm{2}}}}}\left[-\right]\text{(18)}& {\displaystyle }& {\displaystyle }\mathrm{NSE}=\mathrm{1}-\sum _{i=\mathrm{1}}^{n_i}\sum _{j=\mathrm{1}}^J{\displaystyle \frac{{\left[z({\mathit{x}}_i,t_j)-Z^{\ast }({\mathit{x}}_i,t_j)\right]}^{\mathrm{2}}}{{\left[z({\mathit{x}}_i,t_j)-\stackrel{\mathrm{‾}}{z}\right]}^{\mathrm{2}}}}\left[-\right]\text{(19)}& {\displaystyle }& {\displaystyle }\mathrm{B}\mathrm{1}={\displaystyle \frac{\mathrm{1}}{n_{\mathrm{tot}}}}\sum _{i=\mathrm{1}}^{n_i}\sum _{j=\mathrm{1}}^J\left[z({\mathit{x}}_i,t_j)-Z^{*}({\mathit{x}}_i,t_j)\right]\left[\mathrm{mm}\right]\text{(20)}& {\displaystyle }& {\displaystyle }\mathrm{B}\mathrm{2}={\displaystyle \frac{\mathrm{1}}{n_{\mathrm{tot}}}}\sum _{i=\mathrm{1}}^{n_i}{\left[\sum _{j=\mathrm{1}}^J\left(z({\mathit{x}}_i,t_j)-Z^{*}({\mathit{x}}_i,t_j)\right)\right]}^{\mathrm{2}}\left[{\mathrm{mm}}^{\mathrm{2}}\right]\text{(21)}& {\displaystyle }& {\displaystyle }\mathrm{B}\mathrm{3}={\displaystyle \frac{\mathrm{1}}{n_{\mathrm{tot}}}}\sum _{j=\mathrm{1}}^J{\left[\sum _{i=\mathrm{1}}^{n_i}\left(z({\mathit{x}}_i,t_j)-Z^{*}({\mathit{x}}_i,t_j)\right)\right]}^{\mathrm{2}}\left[{\mathrm{mm}}^{\mathrm{2}}\right]\text{(22)}& {\displaystyle }& {\displaystyle }\mathrm{RMSE}=\sqrt{{\displaystyle \frac{\mathrm{1}}{n_{\mathrm{tot}}}}\sum _{i=\mathrm{1}}^{n_i}\sum _{j=\mathrm{1}}^J{\left[z({\mathit{x}}_i,t_j)-Z^{*}({\mathit{x}}_i,t_j)\right]}^{\mathrm{2}}}\left[\mathrm{mm}\right]\end{array}$$

where J is the number of time steps, n_i is the number of observation locations and n_tot is the total number of cross-validated observations. Note that the cross-validation for only one time step (J=1) would yield the following relations: $n_{\mathrm{tot}}\times \mathrm{B}{\mathrm{1}}^{\mathrm{2}}=\mathrm{B}\mathrm{3}$ and B2=RMSE².

3.2.1 Summary results

The overall values of the six objective functions from all 32 226 original observations, along with a separation into dry season (calender months: 4–9) and wet season (calender months: 1–3 and 10–12) are given in Table 1.

Table 1 Summary results from the cross-validation of the estimator Z^∗ for the 12 calender months of the entire year, and split into dry (calender months: 4–9) and wet (calender months: 1–3 and 10–12) season.

Figure 9 Errors of the estimator Z^∗(x₁₂₅) at raingauge “Wilgervier”: evolution of temporal bias B2 over the study period (a) and smoothed distribution of the relative estimation error ϵ_r (b).

The total values of the correlation coefficient ρ, the NSE coefficient, the temporal bias B2 and the RMSE are better for QK-Γ and QK-Wei than for OK and EDK, evoking from a superior performance especially during the wet season when not many of many zero values are present (see Table 1).

Complementary, OK and EDK have superior values for the biases B1 and B3 as a result of the implicit definition as best linear and unbiased estimator. OK (and to some extent EDK) optimize the spatial bias B3 for a given month by adapting their global mean to the observed mean, according to Eq. (1) (Eq. 3). However, this evokes a systematic underestimation in regions of high precipitation and a systematic overestimation in regions of low precipitation. Therefore, a temporal bias B2 accumulates for a location, which consistently experiences extreme precipitation over time. Especially during the wet season, QK outperforms OK and EDK with respect to the temporal bias. The following investigations on raingauge “Wilgervier” exemplary serve as illustration for the evolution of a temporal bias.

3.2.2 Temporal bias at raingauge “Wilgervier”

Raingauge “Wilgervier” (i=125, see Fig. 3) records a relatively high monthly precipitation of 70.1 mm in average compared to the average monthly precipitation of 54.7 mm in the entire domain.

Table 2 Summary results from the cross-validation of the estimation error for the entire year, and split into dry (calender months: 4–9) and wet (calender months: 1–3 and 10–12) season.

The evolution of the temporal bias B2 at raingauge “Wilgervier” is calculated from cross-validation according to Eq. 20 and illustrated in Fig. 9 (left). In addition, the relative estimation error ϵ_r is estimated from the 218 (out of the 264 possible) original observations at “Wilgervier”, being defined as

$$\begin{array}{}\text{(23)}& {\mathit{ϵ}}_r({\mathit{x}}_{\mathrm{125}},t_j)={\displaystyle \frac{Z^{\ast }({\mathit{x}}_{\mathrm{125}},t_j)-z({\mathit{x}}_{\mathrm{125}},t_j)}{z({\mathit{x}}_{\mathrm{125}},t_j)}}.\end{array}$$

The 218 values of ϵ_r(x₁₂₅) are smoothed by a Gaussian kernel with a defined range d_G (=0.35). The distribution of the relative estimation errors should ideally be symmetrical around zero. However, the respective distributions are truncated due the confinement to ${\mathit{ϵ}}_r\ge -\mathrm{1}$ for the variable of monthly precipitation. The smoothed distributions are factually a summary of the estimation errors and are illustrated in Fig. 9 (right).

Figure 10 Evolution of two objective functions for the estimator over the 12 calendar months: Correlation coefficient ρ (a) and RMSE (b). Note that the dry ratio (dashed brown line) is indicated as a percentage on the right axis.

OK displays the highest systematic underestimation over time (Fig. 9 left) and the relative estimation errors have a mode of −20 % (Fig. 9 right). EDK slightly improves the systematic bias of the interpolation, but the relative error distribution still possess a mode of −15 %. QK-Γ and QK-Wei can further improve the systematic underestimation over time and exhibit error distributions with modes of −12 %, and −10 % respectively.

Raingauge “Wilgervier” illustrates that OK and EDK might optimize the spatial bias (Table 1), but they are hampered to minimize the temporal bias in locations of extreme observations. QK, as a spatio-temporal interpolation method, is capable of reducing the temporal bias within regions of relatively high (or low) precipitation, which is potentially important for possible successive water balance considerations.

3.2.3 Cross-validation for different calendar months

The effects of the increased occurrence of zero-value observations on the Pearson correlation coefficient ρ (Eq. 17) and the RMSE (Eq. 22) is exemplary examined next. The respective values are calculated for each calender month from the cross-validation and are illustrated in Fig. 10 along with the dry ratio (Fig. 4).

QK-Γ and QK-Wei display improved values in comparison to OK and EDK for the two selected objective functions from October to March (Fig. 10). However, their performance deteriorates from May to September, when many zero-value observations are present, indicated by a dry ratio of at least 25 % or above. The correlation coefficient ρ plunges in July for both versions of QK (Fig. 10 left) and the respective RMSE shows a similar qualitative behavior (Fig. 10 right).

The performance of QK is considerably influenced by the dry ratio. The presence of many zero values in the data leads to very steep or nearly vertical theoretical cdfs, hampering the allocation of the quantiles to the respective precipitation values.

Figure 11 Histograms for the $F_{Z^{\ast }}$-values of four different interpolation methods.

3.3 Cross-validation of the uncertainty

The estimated error distribution of the estimator Z^∗(x) is described by the associated standard deviation σ_K(x) as a measure of associated uncertainty. The quality of the uncertainty from the cross-validation is assessed by two objective functions: the adapted linear error in the probability space LEPS (Ward and Folland, 1991) and a test on uniformity (Bárdossy and Li, 2008).

LEPS compares the values of the estimator $Z^{\ast }({\mathit{x}}_i,t)$ and the observation z(x_i,t) within the estimated cdf $F_{Z^{\ast }}$ of the error distribution as

$$\begin{array}{}\text{(24)}& \mathrm{LEPS}={\displaystyle \frac{\mathrm{1}}{n_{\mathrm{tot}}}}\cdot \sum _{i=\mathrm{1}}^{n_{\mathrm{tot}}}\left|F_{Z^{\ast }}\left(z\right({\mathit{x}}_i,t\left)\right)-F_{Z^{\ast }}\left(Z^{\ast }\right({\mathit{x}}_i,t\left)\right)\right|,\end{array}$$

LEPS is defined on the interval [0,1]: low values indicate a higher probability for the observation to originate from the estimated probability density distribution and vice versa. The average over the differences of all observations n_tot yields the overall LEPS value.

The test on uniformity verifies the estimated, conditional distribution $F_{Z^{\ast }}$ by calculating its value $F_{Z^{\ast }}\left(z\right({\mathit{x}}_i,t\left)\right)$ for every original observation z(x_i,t). The resulting values (or quantiles) should be uniformly distributed on the interval [0,1] (Bárdossy and Li, 2008). We defined 10 classes of equal width, which should have the same resulting relative frequency. The deviation from uniformity is quantified by the χ²- test variable as the sum of the relative squared differences between uniformity and empirical distribution, ranging from zero (perfect) to nine (improper).

3.3.1 Summary results

The values of the two objective functions from cross-validation of all 32 226 original observations of the entire year, and divided into dry (calender months: 4–9) and wet season (calender months: 1–3 and 10–12) are displayed in Table 2.

The best overall LEPS values are received from the traditional EDK and OK (Table 2). QK-Wei is superior to QK-Γ, but both versions of QK are displaying higher LEPS values than OK or EDK, originating from the dry season when many zero values are present in the data.

Figure 12 Evolution of the two objective functions for the error distribution over 12 calendar months: LEPS (a) and χ² (b). Note that the dry ratio (dashed brown line) is indicated as a percentage on the right axis.

However, the χ²-test variables (Table 2) exhibit a reverse hierarchy among the implemented interpolation methods: QK is superior during the dry season and similar during the wet season. The χ²-test variables should be read in conjunction with the corresponding histograms of the $F_{Z^{\ast }}$ values (Fig. 11). Note that the outer classes of the histograms host all the observations z(x_i,t), which are situated outside the estimated distribution. These classes exhibit the largest deviation from the ideal uniform distribution.

QK-Γ and QK-Wei provide in general a more uniform distribution of $F_{Z^{\ast }}$ than OK and EDK (Table 2), which possess the same value of the χ²-test variable and similar histograms (Table 2 and Fig. 11) due to their implicit affinity in definition.

3.3.2 Cross-validation for different calendar months

The effect of many zero-value observations on the error distribution is investigated by the differentiation into calendar months. The objective functions are recalculated accordingly and illustrated in Fig. 12.

The temporal evolution of the LEPS values for the two versions of QK is influenced by the presence of many zero-value observations. QK-Γ and QK-Wei exhibit LEPS values superior to OK and EDK from September to April, characterized by a dry ratio of less than 26 % (Fig. 12, left). However, the performance of QK deteriorates from May to August when many zero-value observations are present. This dependence explains the overall inferior LEPS values for QK in Table 2. The LEPS values for OK and EDK are hardly influenced by the dry ratio (Fig. 12 left) and show a congruent behavior.

The temporal evolution of the χ² test variable (Fig. 12 right) shows better values for QK-Γ and QK-Wei than for OK or EDK during most calender months. QK maintains a more uniform distribution of the $F_{Z^{\ast }}$ values even for the months with a high dry ratio when OK and EDK deteriorate.

The cross-validation for the uncertainty suggests an improvement by QK under the prerequisite of a low dry ratio within the input data. This improvement is attributed to the wider range of the error distribution and the increased relation between the magnitude of the estimator and the spread of the distribution (see Sect. 3.1).

4 Conclusions

The geostatistical interpolation method of QK addresses the spatial non-stationarity of a variable of interest by its conversion into quantiles and defining distribution parameters. The spatial–temporal description of the variable by QK is a novelty in applied geostatistics and can be regarded as a temporal extension of probability kriging. Therefore, the proposed method could be extended to spatially aggregated variables of streamflows, requiring, however, further investigations. The proposed method accommodates skewed marginal distributions and converts them into an ideal Gaussian distribution prior to interpolation as a major theoretical advantage over the traditional OK or EDK. QK describes an asymmetrical distribution of the random variable Z(x) by the nonlinear estimator Z^∗(x) and the estimation variance ${\mathit{\sigma }}_{\mathrm{K}}^{\mathrm{2}}\left(\mathit{x}\right)$ of the error. QK further establishes a relation between the magnitude of both descriptors.

The variable of monthly precipitation, observed at 226 raingauges over 264 consecutive time steps, serves as input data. We selected the two-parametric Γ distribution and Weibull distribution, because they are defined on the interval [0,∞] and are suitable to describe the variable of monthly precipitation. The selected distributions are fitted to the observations of a specific calendar month, implying an absence of temporal dependence between two sample members (e.g., between the monthly precipitation of December 2002 and December 2003). However, QK does accommodate temporal independence between consecutive observations, unlike existing spatio-temporal kriging methods. In general, other types of distributions, with a higher number of parameters could be selected, especially in case of other variables of interest. Finally, we used elevation as external drift, both for the interpolation of the parameters within QK as well as for the reference EDK.

The cross-validation of the estimator revealed an improvement for most of the selected objective functions. In particular, QK addresses the temporal bias, which remains unattended by the traditional geostatistical methods, which only optimize the mean spatial bias. In case of the estimator, QK-Γ performs slightly better than QK-Wei for most of the selected objective functions. The cross-validation of the associated uncertainty shows an improvement by QK in the description of the distribution of the estimation errors in comparison to the traditional geostatistical interpolation methods. However, its performance depends on the percentage of zero values in the input data and declines when many zero values are present. In general, QK-Wei shows a superior estimation of the associated uncertainty than QK-Γ.