Conditional simulation of surface rainfall fields using modified phase annealing

The accuracy of quantitative precipitation estimation (QPE) over a given region and period is of vital importance across multiple domains and disciplines. However, due to the intricate temporospatial variability and the intermittent nature of precipitation, it is challenging to obtain QPE with adequate accuracy. This paper aims to simulate rainfall fields while honoring both the local constraints imposed by the point-wise rain gauge observations and the global constraints imposed by the field measurements obtained from weather radar. The conditional simulation method employed in this study is modified phase annealing (PA), which is practically an evolution from the traditional simulated annealing (SA). Yet unlike SA, which implements perturbations in the spatial field, PA implements perturbations in Fourier space, making it superior to SA in many respects. PA is developed in two ways. First, taking advantage of the global characteristic of PA, the method is only used to deal with global constraints, and the local ones are handed over to residual kriging. Second, except for the system temperature, the number of perturbed phases is also annealed during the simulation process, making the influence of the perturbation more global at initial phases and decreasing the global impact of the perturbation gradually as the number of perturbed phases decreases. The proposed method is used to simulate the rainfall field for a 30 min event using different scenarios: with and without integrating the uncertainty of the radar-indicated rainfall pattern and with different objective functions.

1 Introduction

Quantitative precipitation estimation (QPE) over a given region and period is of vital importance across multiple domains and disciplines. Yet the intricate temporospatial variability, together with the intermittent nature of precipitation in both space and time, has hampered the accuracy of QPE (Kumar and Foufoula-Georgiou, 1994; Emmanuel et al., 2012; Cristiano et al., 2017).

The point-wise observations of precipitation measured by rain gauges are accurate but only available at limited locations. Meanwhile, precipitation-related measurements produced by meteorological radars have become standard outputs of weather offices in many places in the world. However, the problem with radar-based QPE is the nonguaranteed accuracy, which could be impaired by various sources of errors, such as static/dynamic clutter, signal attenuation, anomalous propagation of the radar beam, uncertainty in the Z–R relationship, etc. (see Doviak and Zrnić, 1993; Collier, 1999; Fabry, 2015, for details). Despite the various sources of error, weather radar has been widely acknowledged as a valid indicator of precipitation patterns (e.g., Mendez Antonio et al., 2009; Fabry, 2015). Considering the pros and cons of the two most usual sources of precipitation information, the QPE obtained by merging the point-wise rain gauge observations and the radar-indicated precipitation pattern has become a research problem in both meteorology and hydrology (Hasan et al., 2016; Yan and Bárdossy, 2019).

In the context of merging radar and rain gauge data, we consider two types of constraints: the local constraints imposed by the point-wise rain gauge observations and the global constraints imposed by the field measurement from weather radar. This paper focuses on simulating surface rainfall fields conditioned on the two types of constraints. There exists a variety of geostatistical methods aiming to simulate conditional Gaussian fields with a given covariance function, such as turning bands simulation, LU-decomposition-based methods, sequential Gaussian simulation, etc. (see Deutsch and Journel, 1998; Chilès and Delfiner, 2012; Lantuéjoul, 2013, for details). The common goal of these methods is to ensure that the simulated realizations comply with the additional information available (Lauzon and Marcotte, 2019). The additional information could be observed values of the simulated targets, measurements that are related linearly or nonlinearly to the simulated targets, third- or higher-order statistics (Guthke and Bárdossy, 2017; Bárdossy and Hörning, 2017), etc.

The conditional simulation method used in this work is phase annealing (PA). It was first proposed in Hörning and Bárdossy (2018) and is essentially a product of the general-purpose metaheuristics method of simulated annealing (SA) (Kirkpatrick et al., 1983; Geman and Geman, 1984; Deutsch, 1992; Deutsch and Journel, 1994). It utilizes the sophisticated optimization scheme of SA in the search for the global optimum. Yet compared to  SA, the distinction or evolvement of PA lies in the fact that the perturbation, or swapping in the nomenclature of SA, is implemented in Fourier space. Or, to be more exact, the perturbation is implemented on the phase component of the Fourier transform, while the power spectrum is preserved such that the spatial covariance is invariant at all iterations according to the well-known Wiener–Khinchin theorem (Wiener, 1930; Khintchine, 1934). Compared to SA, PA alleviates the singularity problem, namely the undesired discontinuities or poor embedding of the conditional points within the neighborhood (Hörning and Bárdossy, 2018), and in general, PA has a much higher convergence rate compared to SA.

A remarkable feature of PA is that it is a global method: any perturbation imposed on the phase component is reflected in the entire field. Yet admittedly, if the perturbation is implemented at lower frequencies where the corresponding wave lengths are relatively large, the impact is more global and vice versa. The global characteristic of PA makes it a valuable methodology for global constraints. However, PA is found to be insufficient in the treatment of local constraints (Hörning and Bárdossy, 2018). Note that by local constraints, we refer primarily to point equality constraints when the total number of constraints is far less than that of the grid points. In the algorithm of PA, the local constraints are normally ensured by inserting a component measuring the dissimilarity between the simulated and the target values. On the other hand, one could argue that if the information on the measurement error is explicitly known, then this piece of information could be considered in the simulation, and the inability of PA to handle the local constraints can be utilized in turn. However, in the general case, the local constraints can only be approximated by PA.

Respecting the fact that the specialty of PA is the treatment of the global constraints, we separate the global from the local. In particular, PA is only used to handle the global constraints, and local ones are handled separately by residual kriging at each iteration. As an extension of PA, except for annealing the system temperature, the number of perturbed phases is annealed in parallel to render the algorithm to work more globally at initial phases of the simulation. The global impact of the perturbation is weakened as the number of perturbed phases decreases.

This paper is divided into six sections. After the general introduction, the methodology of PA is introduced in Sect. 2, including three stages: presimulation, simulation, and postsimulation. Section 3 provides two options to integrate the uncertainty of the radar-indicated rainfall pattern into the simulation. Section 4 introduces the study domain and the two types of data used in this study. In Sect. 5, the proposed algorithm is used to simulate the rainfall field for a 30 min event, where different simulation scenarios are applied. Section 6 ends this paper with conclusions.

2 Methodology

Figure 1 summarizes the procedure of simulating surface rainfall fields using the algorithm of PA, including three stages: presimulation (PreSim), simulation (Sim), and postsimulation (PostSim). Each stage and the corresponding substages are described in the following subsections in the same sequence as shown in the flowchart.

Figure 1 Flowchart of the procedure to simulate surface rainfall fields using the algorithm of PA.

2.1 Presimulation

2.1.1 Distribution function of surface rainfall

PA is embedded in Gaussian space. The direct output of the PA algorithm is a spatial field whose marginal distribution function is standard normal; hence the distribution function of surface rainfall is essential to transform the simulated Gaussian fields into rainfall fields of interest.

The scenario is as follows: based on K rain gauge observations and a regular grid of radar quantiles (representing the spatial pattern of the rainfall field), the distribution function of surface rainfall is generated. The procedure was first introduced in Yan and Bárdossy (2019), and we specify the modified version here.

• a.

  For all rain gauge observations r_k, the collocated quantiles in the radar quantile map U are determined and denoted as u_k. The two datasets are then sorted in ascending order, i.e., $r_{\mathrm{1}}\le \mathrm{\cdots }\le r_K$ and $u_{\mathrm{1}}\le \mathrm{\cdots }\le u_K$.

• b.

  Set the spatial intermittency u₀, i.e., the quantile corresponding to zero precipitation. u₀ is estimated as the ratio of the number of zero-valued pixels to the total number of pixels in the domain of interest. If the small quantiles of the field are properly sampled by the rain gauges, one could use the smallest sampled quantile u₁ to estimate u₀:

  $$\begin{array}{}\text{(1)}& u_{\mathrm{0}}=\left\{\begin{array}{ll}{u}_{\mathrm{1}}& \mathrm{if}{r}_{\mathrm{1}}=\mathrm{0}\\ u_{\mathrm{1}}/\mathrm{2}& \mathrm{otherwise},\end{array}\right.\end{array}$$

  where r₁ is the smallest gauge observation corresponding to u₁. As an alternative, if the small quantiles of the field are not properly sampled by the rain gauges, we propose to estimate u₀ from the original radar displays in dBZ, namely a series of instantaneous maps of radar reflectivity within the relevant accumulation time. A zero-valued pixel is defined when the maximum of the collocated pixel values of these maps is smaller than a given threshold (say 20 dBZ).

• c.

  Let G(r) denote the distribution function of surface rainfall to be estimated. A linear interpolation is applied for rainfall values of less than the maximum recorded gauge observations, i.e., r≤r_K:

  $$\begin{array}{}\text{(2)}& G\left(r\right)={\displaystyle \frac{u_k-u_{k-\mathrm{1}}}{r_k-r_{k-\mathrm{1}}}}\left(r-r_{k-\mathrm{1}}\right)+u_{k-\mathrm{1}},\end{array}$$

  with r_k−1, r_k being the two nearest neighbors of r ($r_{k-\mathrm{1}}\le r\le r_k$) and $u_{k-\mathrm{1}},u_k$ being the quantiles corresponding to r_k−1, r_k, respectively.

• d.

  Extrapolate rainfall values r>r_K. A modification is made here: the minimum of exponential and linear extrapolation is used as the result, as expressed in Eq. (3). We have learned from practice that the exponential extrapolation tends to overestimate the rainfall extremes, so a linear component is used to restrict the extrapolation result.

  $$\begin{array}{}\text{(3)}& G\left(r\right)=\min \left(\mathrm{1}-e^{-\mathit{\lambda }r},{\displaystyle \frac{u_K-u_{K-\mathrm{1}}}{r_K-r_{K-\mathrm{1}}}}\left(r-r_K\right)+u_K\right),\end{array}$$

  where the only parameter λ of the exponential distribution is determined from the last pair (r_K, u_K):

  $$\begin{array}{}\text{(4)}& \mathit{\lambda }=-\ln \left(\mathrm{1}-u_K\right)/r_K.\end{array}$$

It is recommended that the methodology described above is used to estimate the distribution function of accumulated rainfall, not rain intensity, because it is built on the assumption that rain gauge observations can represent the ground truth. Yet in general, rain gauges are considered to be poor in capturing the instantaneous rain intensity, while the measurement error diminishes rapidly as the integration time increases (Fabry, 2015).

As we only use the radar-indicated spatial ranks, i.e., the scaled radar map following a uniform distribution in [0, 1], there is no requirement for the accuracy of the Z–R relation given the monotonic relationship of the two quantities. One could use a Z–R relation at hand to transform radar reflectivity into rain intensity and then make the accumulation. Nevertheless, a normal quality control for radar data, as described in Sect. 4, is necessary to maintain the accuracy of radar-indicated spatial ranks.

2.1.2 Marginal conversion to Gaussian

As PA is embedded in Gaussian space, all the constraints, including the point equality targets r_k (mm) imposed by rain gauges and the quantile map U (–) imposed by weather radar, are converted to a standard normal marginal distribution by

$$\begin{array}{}\text{(5)}& {\displaystyle }& {\displaystyle }{z}_k={\mathrm{\Phi }}^{-\mathrm{1}}\left(G\left(r_k\right)\right)\mathrm{for}k=\mathrm{1},\mathrm{\cdots },K\text{(6)}& {\displaystyle }& {\displaystyle }{Z}^{*}={\mathrm{\Phi }}^{-\mathrm{1}}\left(U\right),\end{array}$$

where Φ⁻¹ is the inverse of the cumulative standard normal distribution function and G is the distribution function of surface rainfall obtained according to the procedure described in Sect. 2.1.1.

2.1.3 Objective function

We impose two kinds of constraints: local and global. As has been explained in the introduction, PA is a powerful method for handling global constraints. In order not to interfere with the logic behind PA and to fulfill the point equality constraints (local constraints) exactly, residual kriging is implemented at each iteration to obtain the observed values at rain gauge locations. Thus the objective function only needs to measure the fulfillment of the global constraints.

We impose two global constraints: field pattern and directional asymmetry. Note that both constraints are evaluated from the simulated Gaussian field and are compared with the radar-based Gaussian field Z^*, as obtained in Eq. (6). We term Z^* the reference field hereafter.

The first global constraint, field pattern, requires the simulated field to be similar to the reference field. The similarity of the two fields is quantified by the Pearson correlation coefficient:

$$\begin{array}{}\text{(7)}& {\mathit{\rho }}_{Z,Z^{*}}={\displaystyle \frac{\mathrm{cov}\left(Z,Z^{*}\right)}{{\mathit{\sigma }}_Z{\mathit{\sigma }}_Z^{*}}}.\end{array}$$

In the ideal case, ${\mathit{\rho }}_{Z,Z^{*}}$ equals 1, and we use the difference ($\mathrm{1}-{\mathit{\rho }}_{Z,Z^{*}}$) to measure the distance from the ideal.

The second global constraint is directional asymmetry, as given by

$$\begin{array}{}\text{(8)}& \left[A\right(h){]}_Z={\displaystyle \frac{\mathrm{1}}{N\left(h\right)}}\sum _{x_i-x_j\approx h}{\left(\mathrm{\Phi }\left(Z\left(x_i\right)\right)-\mathrm{\Phi }\left(Z\left(x_j\right)\right)\right)}^{\mathrm{3}},\end{array}$$

where [A(h)]_Z (abbreviated as A_Z hereafter) is the asymmetry function evaluated from the simulated Gaussian field Z, N(h) is the number of pairs fulfilling $x_i-x_j\approx h$, and Φ is the cumulative standard normal distribution function. Directional asymmetry was first introduced in Bárdossy and Hörning (2017) and Hörning and Bárdossy (2018). It is a third-order statistic, and the physical phenomenon revealed by this statistic could be significant for advection-dominant processes, such as storms. We use the entire field to compute the directional asymmetry function and compare the simulated one with the reference asymmetry function, i.e., the directional asymmetry function evaluated from the reference field, abbreviated as $A_{Z^{*}}$. There exist multiple choices to define the distance of the two asymmetry functions, and we have used the L^∞ norm $\Vert A_Z-A_{Z^{*}}{\Vert }_{L^{\mathrm{\infty }}}$.

Different schemes could be used to combine the two components, linear or nonlinear, and we have chosen the maximum of the two components. Finally, the objective function we have used to quantify the fulfillment of the two global constraints is expressed as

$$\begin{array}{}\text{(9)}& \mathcal{O}\left(Z\right)=\max \left(\left(\mathrm{1}-{\mathit{\rho }}_{Z,Z^{*}}\right),\left(w/A_{\mathrm{scal}}\right)\cdot \Vert A_Z-A_{Z^{*}}{\Vert }_{L^{\mathrm{\infty }}}\right),\end{array}$$

where w is the relative weight of the component directional asymmetry, and A_scal, as expressed in Eq. (10), is the scaling factor that scales the L^∞ norm of the difference of the two asymmetry functions between 0 and 1.

$$\begin{array}{}\text{(10)}& A_{\mathrm{scal}}=\Vert A_{Z^{*}}{\Vert }_{L^{\mathrm{\infty }}}\end{array}$$

It is worth mentioning that the proposed methodology is flexible, and any global statistic could be used to constrain the simulated fields. One could even combine several statistics in the objective function.

2.2 Simulation

2.2.1 Stopping criteria and cooling schedule

There are many choices of stopping criteria for an optimization algorithm, such as (a) the total number of iterations, (b) the predefined value of the objective function, (c) the rate of decrease of the objective function, (d) the number of continuous iterations without improvement, (e) the predefined threshold of the initial objective function value, and so forth. One could use one criterion or combine several.

As for the cooling schedule, the system temperature of PA decreases according to the cooling schedule as the optimization process develops: the lower the system temperature, the less likely a bad perturbation is accepted. A bad perturbation occurs when the perturbation does not decrease the objective function value. A reasonable cooling schedule is capable of preventing the optimization from being trapped prematurely at a local minimum. Yet one should be aware that there is always a compromise between the statistical guarantee of the convergence and the computational cost: the slower the temperature decreases, the higher the probability of the convergence. However, cooling slowly also means more iterations and therefore higher computational costs.

Comparative studies of the performance of SA using the most important cooling schedules, i.e., multiplicative monotonic, additive monotonic, and nonmonotonic adaptive cooling, have been made by, e.g., Nourani and Andresen (1998), Martín and Sierra (2020), and others. The results show that annealing works properly when the cooling curve has a moderate slope at the initial and central stages of the process and tends to have a softer slope at the final stage. Many cooling schedules satisfy these conditions. As PA utilizes the optimization strategy of SA, the rules also apply for PA. Our choice of the cooling schedule is the multiplicative monotonic one.

Specifically, we search for two parameters – the initial temperature T₀ and the final temperature T_min – by decreasing the temperature exponentially and discretely. At each fixed temperature, N perturbations, say 1000, are implemented, and the corresponding acceptance rate and improvement rate are computed. Acceptance means the perturbation is accepted by the system, and the system state is updated irrespective of whether the perturbation brings improvement to the system, while improvement means the perturbation does decrease the objective function value. In short, a perturbation that brings improvement must be an accepted perturbation, yet an accepted perturbation does not necessarily bring improvement to the system. We refer to the experiment as a temperature cycle, where N perturbations are implemented, and the two mentioned rates are computed at a fixed temperature.

T₀ is first set by starting a temperature cycle with an initial guess of T₀. If the acceptance rate is lower than the predefined limit, say 98 %, then the temperature is increased and vice versa. The goal is to find a T₀ with a relatively high acceptance rate. Then, T_min is set by decreasing T₀ exponentially until the predefined stopping criterion is met. It is clear in our case that if more iterations are implemented, better realizations (realizations with lower objective function values) could be produced, yet it is again a compromise between a satisfactory destination and the computational cost.

We note that from the determination of T₀, the total number of temperature cycles m (to anneal T₀ to T_min) is recorded. Thus the total number of perturbations is estimated as L=mN. A continuous cooling schedule is then computed from T₀, T_min, and L. The temperature at Iteration l is computed as

$$\begin{array}{}\text{(11)}& T_l=T_{\mathrm{0}}\cdot {\mathit{\alpha }}_T^l\mathrm{for}l=\mathrm{0},\mathrm{\cdots },L-\mathrm{1},\end{array}$$

where α_T is the attenuation factor of the system temperature, given as

$$\begin{array}{}\text{(12)}& {\mathit{\alpha }}_T=\exp \left(\ln \left({\displaystyle \frac{T_{\min }}{T_{\mathrm{0}}}}\right)/(L-\mathrm{1})\right).\end{array}$$

In parallel with the temperature annealing we also anneal the number of phases being perturbed in Fourier space, starting with a relatively large number N₀ (say 5 % to 20 % of the total number of valid Fourier phases) and decreasing all the way down to 1. The number of phases to be perturbed at Iteration l is computed by

$$\begin{array}{}\text{(13)}& N_l=N_{\mathrm{0}}\cdot {\mathit{\alpha }}_N^l\mathrm{for}l=\mathrm{0},\mathrm{\cdots },L-\mathrm{1},\end{array}$$

where α_N is the attenuation factor with respect to the number of perturbed Fourier phases. It is computed as

$$\begin{array}{}\text{(14)}& {\mathit{\alpha }}_N=\exp \left(-\ln \left(N_{\mathrm{0}}\right)/(L-\mathrm{1})\right).\end{array}$$

The logic behind annealing the number of perturbed phases is to render the perturbation to have a more global impact initially when the distance from the destination is relatively large and to decrease the impact of the perturbation gradually as the optimization process develops.

2.2.2 Starter generation

PA requires a starting Gaussian random field with the prescribed spatial covariance, abbreviated as starter. Various methods could be used to generate such a random field, e.g., fast Fourier transformation for regular grids (Wood and Chan, 1994; Wood, 1995; Ravalec et al., 2000), turning bands simulation (Journel, 1974), or the Cholesky transformation of the covariance matrix. Considering the fact that the gauge observations used in this work cannot provide enough information to derive the spatial covariance, the radar-based Gaussian field, i.e., the reference field Z^*, is used instead to derive the covariance.

It should be noted that if a fast Fourier transformation (FFT) is used to generate the starter, the inherent periodic property of the FFT should be treated with care. Specifically, the simulation should be embedded in a larger domain. The original domain is enlarged in all directions by a finite range, i.e., the range bringing the covariance function from the maximum to approaching zero. If covariance models with asymptotic ranges (e.g., exponential, Matérn, and Gaussian covariances) are employed, the extension in domain size could be significant (Chilès and Delfiner, 2012).

2.2.3 Main PA cycle

The starting point of PA is a Gaussian random field Z_K with the prescribed spatial covariance as described in Sect. 2.2.2. A little modification in our case is that the values at the observational locations are fixed by residual kriging, as indicated by the subscript K. In particular, the procedure of residual kriging is explained in the following (Step 6). The procedure of the PA algorithm applied in this paper is specified as follows:

1. The discrete Fourier transform (DFT) of Z_K is computed as

   $$\begin{array}{}\text{(15)}& {\mathcal{Z}}_K=\mathcal{F}\mathit{\left\{}{Z}_K\mathit{\right\}}.\end{array}$$

2. The system temperature T_l and the number of perturbed phases N_l at Iteration l are computed by Eqs. (11) and (13).

3. N_l vectors (u_n, v_n) are generated, where n=0, ⋯, N_l−1. u_n, v_n are randomly drawn from the two discrete uniform distributions [1, ⋯, U^′] and [1, ⋯, V^′], respectively, where U^′, V^′ are the highest frequencies considered for perturbation in the x and y directions, respectively. Note that zero frequencies in both directions should be excluded from the perturbation. One could use the entire frequency range or select a subarea to impose the perturbation.

4. N_l phases θ_n are randomly drawn from the uniform distribution [−π, π).

5. The Fourier coefficients at the selected locations (u_n, v_n) and the corresponding symmetrical locations are expressed as

   $$\begin{array}{}\text{(16)}& {\displaystyle }& {\displaystyle }{\mathcal{Z}}_K\left[u_n,v_n\right]=a_n+j{b}_n\text{(17)}& {\displaystyle }& {\displaystyle }{\mathcal{Z}}_K\left[U-u_n,V-v_n\right]=a_n-j{b}_n,\end{array}$$

   where U, V are the number of grid points in x and y directions, respectively, and are also the highest frequencies in both directions. These coefficients are then updated in terms of the Fourier phases using θ_n, while the corresponding amplitudes remain unchanged:

   $$\begin{array}{}\text{(18)}& {\displaystyle }{\mathcal{Z}}_K\left[u_n,v_n\right]\leftarrow \sqrt{a_n^{\mathrm{2}}+b_n^{\mathrm{2}}}\cdot \left(\cos \left({\mathit{\theta }}_n\right)+j\sin \left(\mathit{\theta }\right)\right)\text{(19)}& {\displaystyle }\begin{array}{rl}{\mathcal{Z}}_K\left[U-u_n,V-v_n\right]& \leftarrow \sqrt{a_n^{\mathrm{2}}+b_n^{\mathrm{2}}}\\ & \cdot \left(\cos \left({\mathit{\theta }}_n\right)-j\sin \left({\mathit{\theta }}_n\right)\right).\end{array}\end{array}$$

   The perturbed Fourier transform is denoted as $\stackrel{\mathrm{̃}}{\mathcal{Z}}$. The corresponding perturbed spatial random field is obtained by the inverse DFT:

   $$\begin{array}{}\text{(20)}& \stackrel{\mathrm{̃}}{Z}={\mathcal{F}}^{-\mathrm{1}}\mathit{\left\{}\stackrel{\mathrm{̃}}{\mathcal{Z}}\mathit{\right\}}.\end{array}$$

6. Due to the perturbation, $\stackrel{\mathrm{̃}}{Z}$ no longer satisfies the point equality constraints (note the removal of the subscript K). Thus kriging is applied for the residuals $z_k-\stackrel{\mathrm{̃}}{Z}\left(x_k\right)$, where z_k represents the point equality constraints defined in Eq. (5). The results of kriging r are superimposed on $\stackrel{\mathrm{̃}}{Z}$:

   $$\begin{array}{}\text{(21)}& {\stackrel{\mathrm{̃}}{Z}}_K=\stackrel{\mathrm{̃}}{Z}+r.\end{array}$$

   As kriging is a geostatistical method that depends only on the configuration of the data points, the weight matrix of individual grid points does not change with iterations. One needs to compute the weight matrices for all the grid points only once. Thus residual kriging is cheap to use and causes almost no additional computational cost.

7. ${\stackrel{\mathrm{̃}}{Z}}_K$ is then subject to the objective function defined in Eq. (9). If $\mathcal{O}\left({\stackrel{\mathrm{̃}}{Z}}_K\right)<\mathcal{O}\left(Z_K\right)$, the perturbation is accepted. Otherwise, the perturbation is accepted with the probability

   $$\begin{array}{}\text{(22)}& P=\exp \left({\displaystyle \frac{\mathcal{O}\left(Z_K\right)-\mathcal{O}\left({\stackrel{\mathrm{̃}}{Z}}_K\right)}{T_l}}\right).\end{array}$$

8. If the perturbation is accepted, the system state is updated; namely Z_K, 𝒵_K and 𝒪(Z_K) are replaced by ${\stackrel{\mathrm{̃}}{Z}}_K$, ${\stackrel{\mathrm{̃}}{\mathcal{Z}}}_K$, and $\mathcal{O}\left({\stackrel{\mathrm{̃}}{Z}}_K\right)$, respectively.

9. If the stopping criterion is met, the optimization is stopped and Z_K is a realization in Gaussian space, satisfying the predefined optimization criterion. Otherwise go to Step 2.

2.3 Postsimulation

The Gaussian field Z_K is transformed into a rainfall field by

$$\begin{array}{}\text{(23)}& R=G^{-\mathrm{1}}\left(\mathrm{\Phi }\left(Z_K\right)\right),\end{array}$$

where Φ is the cumulative standard normal distribution function, and G⁻¹ is the inverse of the distribution function of surface rainfall, as obtained in Sect. 2.1.1. The resultant R is a realization of the surface rainfall field.

3 Uncertainty of radar-indicated rainfall pattern

The two datasets – gauge observations r_k and the collocated radar quantiles u_k – are supposed to have the same rankings; namely large gauge observations should correspond to large quantiles and vice versa. Yet it is not always the case, which reveals the conflict of radar and gauge data. To quantify the conflict, we employ the (Spearman's) rank correlation ρ^S that measures the monotonic relationship of two variables: the closer ρ^S is to 1, the more accordance (or less conflict) of radar and gauge data.

The conflict of the two could be partially explained by the fact that weather radar measures at some distance above the ground (a few hundred to more than a thousand meters aloft). It is therefore reasonable to suspect the correctness of comparing the ground-based rain gauge observations with the collocated radar data by assuming the vertical descent of the hydrometeors. In fact, hydrometeors are very likely to be laterally advected during their descent by wind, which occurs quite frequently concurrently with precipitation. To take the possible wind-induced displacement into account, the procedure described in Yan and Bárdossy (2019) is adopted. Here we recall the procedure briefly in three steps.

1. A rank correlation matrix ${\mathit{\rho }}_{ij}^{\mathrm{S}}$ is generated by first shifting the original radar quantile map U with all the vectors defined by a regular grid h_ij, i.e., locations of all the grid cells in h_ij, and then calculating the rank correlation between gauge observations and the collocated radar quantiles in the shifted map. Note that the radar quantile map shifted by vector h_ij (a certain grid cell in h_ij) is denoted as U_ij.

   Both h_ij and ${\mathit{\rho }}_{ij}^{\mathrm{S}}$ have the same spatial resolution as U. The center of h_ij is (0, 0), with the corresponding entry in ${\mathit{\rho }}_{ij}^{\mathrm{S}}$ denoted as ${\mathit{\rho }}_{\mathrm{00}}^{\mathrm{S}}$, meaning that zero shift is imposed on U. The farther a grid cell is from the center, the larger the shift imposed on U is. One should limit the number of grid cells in h_ij depending on the radar measurement height.

2. A probability matrix P_ij is generated from the rank correlation matrix ${\mathit{\rho }}_{ij}^{\mathrm{S}}$:

   $$\begin{array}{}\text{(24)}& {\mathrm{P}}_{ij}=\left\{\begin{array}{ll}\mathrm{0}& \mathrm{for}{\mathit{\rho }}_{ij}^{\mathrm{S}}\le {\mathit{\rho }}_{\mathrm{00}}^{\mathrm{S}}\\ H\left({\mathit{\rho }}_{ij}^{\mathrm{S}}\right)& \mathrm{otherwise},\end{array}\right.\end{array}$$

   where H is a monotonic function, such as power function, exponential function, or logarithmic function. The first derivative (dH∕dx) matters as a large dH∕dx means that more weights are assigned to those displacement vectors (or shifted fields) that produce a higher-rank correlation. Our choice of function H is simply H(x)=x².

3. The probability matrix P_ij is scaled to ensure that the sum of all entries equals 1:

   $$\begin{array}{}\text{(25)}& {\mathbf{P}}_{ij}={\displaystyle \frac{P_{ij}}{\sum _i\sum _j{P}_{ij}}}.\end{array}$$

   P_ij quantifies the uncertainty of the radar-indicated rainfall pattern by giving the probability of individual shifted quantile fields in analogy with a probability mass function. The logic behind it is that only those displacement vectors that increase the accordance of radar and gauge data are accepted and given strong probability.

The method described above is built on the assumption that the entire field is displaced by a single vector. The extent of the domain in this work is 19 km × 19 km (see Sect. 4 for details), and the selection of the extent is limited by the gauge data availability. With this spatial extent, we could quite often find dozens of displacement vectors that bring up the concordance of radar and gauge data. If gauge data accessibility is not the limiting factor, one could enlarge the extent (say to 40 km × 40 km), provided that the resultant probability matrix has some pattern (not a random field).

There are two options to integrate the uncertainty of the radar-indicated rainfall pattern, i.e., the information carried by the probability matrix P_ij, into the simulation.

Option 1: the expectation of the shifted quantile fields is computed using the probability matrix, and then the corresponding Gaussian field is computed, as given in Eq. (26). The result is used as the reference field when applying the PA algorithm.

$$\begin{array}{}\text{(26)}& Z^{*}={\mathrm{\Phi }}^{-\mathrm{1}}\left(\sum _i\sum _j{P}_{ij}\cdot U_{ij}\right)\end{array}$$

Option 2: those (marginally converted) shifted quantile fields Φ⁻¹(U_ij) bearing a positive probability are considered as reference fields, and the algorithm described in Sect. 2 is applied independently for all these reference fields. The results are multiple realizations of the surface rainfall field with different references, denoted as R_ij. The estimate of the surface rainfall field (with the uncertainty of the radar-indicated rainfall pattern integrated) is obtained as the expectation of these realizations:

$$\begin{array}{}\text{(27)}& \stackrel{\mathrm{‾}}{R}=\sum _i\sum _j{P}_{ij}\cdot R_{ij}.\end{array}$$

Both options are capable of producing estimates of the surface rainfall field with the uncertainty of the radar-indicated rainfall pattern integrated, yet the results are distinct, as presented in Sect. 5.

4 Study domain and data

The study domain is located in Baden-Württemberg, in the southwest of Germany. As shown in Fig. 2, it is a square domain with a side length of 19 km. The domain is discretized to a 39×39 grid with a spatial resolution of 500 m × 500 m. A gauge network consisting of 12 pluviometers is available within the domain, as denoted by the red dots in Fig. 2.

Figure 2 The elevation map of Baden-Württemberg with the study domain marked by the red square, the rain gauges marked by the red dots, and the site of Radar Türkheim marked by the red star.

The radar data used in this study are the raw data with 5 min temporal resolution, measured by Radar Türkheim, a C-band radar operated by the German Weather Service (DWD). Radar Türkheim is located about 45 km from the domain center, as denoted by the red star in Fig. 2. The raw data are subject to a processing chain consisting of (1) clutter removal (Gabella and Notarpietro, 2002), (2) attenuation correlation (Krämer and Verworn, 2008; Jacobi and Heistermann, 2016), (3) reprojection from polar to Cartesian coordinates, and (4) clipping of the square data for the study domain. All these quality control steps are operated in the wradlib environment, an open-source library for weather radar data processing (Heistermann et al., 2013).

Figure 3 (a) Original radar quantile map (–) for the 30 min event (4 May 2017, 13:20–13:50 LT), with the rain gauge observations (mm) labeled in yellow and the collocated radar quantiles labeled in cyan. (b) Distribution functions of surface rainfall. Black: original radar map without displacement; gray: radar maps with displacement; red: the weighted average of the radar maps with displacement.

5 Application

We selected a 30 min event to apply the algorithm of PA. The event was selected not only due to the relatively prominent rain intensity reflected by the rain gauge data, but more importantly because the event was properly captured by a few rain gauges unevenly distributed in the domain of interest, as shown by the red dots in Fig. 3a. From the figure, it can be clearly seen that the sampled quantiles (indicated by the text in cyan) cover the entire range [0, 1]: not only the small or large ones, but the sampled quantiles as well are more or less evenly distributed. One might have noticed that the smallest sampled quantile is 0.53. Yet in this case, u₀=0.26 (the spatial intermittency) and the quantile value 0.53 correspond to the rainfall value 0.28 mm after the re-ordering.

The conflict of radar and gauge data is obviously reflected in Fig. 3 (left), for example, the collocation of 4.20 mm rainfall with the quantile 0.99 in contrast with the collocation of 8.24 mm rainfall with the quantile 0.94. The rank correlation of the gauge observations and the collocated radar quantiles is 0.601, as labeled in the title of the figure. The distribution function of the surface rainfall field based on the two original datasets is shown by the black line in Fig. 3b.

However, using the algorithm described in Sect. 3 to evaluate the uncertainty of the radar-indicated rainfall pattern, one can obtain multiple distribution functions of surface rainfall, as shown by the gray lines in Fig. 3b. Note that only the distribution functions associated with those shifted fields possessing a positive probability (as indicated by the probability matrix) are shown in the figure.

With these distribution functions, one can transform the corresponding shifted quantile fields into rainfall fields, and with the probability matrix, the expectation of these rainfall fields can be computed, as shown in Fig. 4b. Also with the probability matrix, one can compute the expectation of the shifted quantile fields, as shown in Fig. 4a. From the expected quantile map, a tangible improvement in terms of the concordance of radar and gauge data is observed: the rank correlation increases from 0.601 to 0.769, as labeled in the title of Figs. 3 and 4.

Figure 4 (a) The expected quantile map (–), with the rain gauge observations (mm) labeled in yellow and the collocated radar quantiles labeled in cyan. (b) The expected rainfall field (mm).

It is noteworthy that the data configuration used in this work is not good enough to maximize the performance of the proposed methodology as the distribution of the rain gauges is relatively centered. Thus we have to select events whose storm center is relatively centered according to the radar map. To enlarge the applicability, it is recommended that rain gauges are uniformly distributed in the domain of interest. In addition to the distribution of the rain gauges, as shown in the following, the effect of conditioning on gauge observations is local (in consistency with the terminology local constraints); hence, with a denser network of rain gauges, a good performance of the methodology is more likely guaranteed.

We open up two simulation sessions, depending on the different objective functions used when applying the PA algorithm. In Session I, the objective function contains only the component field pattern $\mathcal{O}\left(Z\right)=\mathrm{1}-{\mathit{\rho }}_{Z,Z^{*}}$. In Session II, the component directional asymmetry comes into play, and the objective function expressed in Eq. (9) is employed, where the relative weight of the component directional asymmetry w equals 0.5. Technically, for both simulation sessions, the optimization process stops when the objective function falls below 0.05.

5.1 Simulation Session I

We present an evolution in terms of the simulation strategy where the algorithm of PA is applied differently in three stages:

1. using the original quantile map as the reference

2. using the expected quantile map as the reference, i.e., integrating the uncertainty of the radar-indicated rainfall pattern via Option 1 in Sect. 3

3. simulating independently using those shifted quantile maps with a positive probability as the reference and computing the expectation, i.e., integrating the uncertainty of the radar-indicated rainfall pattern via Option 2 in Sect. 3.

Figure 5 Stage 1 results: mean and standard deviation of 100 realizations (mm), obtained using the original quantile map as the reference.

Stage 1 (simulating rainfall fields using the original quantile map as the reference) means the uncertainty of the radar-indicated rainfall pattern is not integrated. Figure 5 shows the mean of 100 such realizations on the left and the corresponding standard deviation map on the right. The standard deviation map reflects the variability of different realizations, and the pixel values, collocated with the rain gauges, are zeros, as marked by the white points in Fig. 5 (right), which shows the fulfillment of the point equality constraints and the local effect of conditioning to the gauge observations. And this local effect is reflected in all the standard deviation maps shown in the following.

Stage 2 (simulating rainfall fields using the expected quantile map as the reference) integrates the uncertainty of the radar-indicated rainfall pattern via Option 1, as described in Sect. 3. Similarly, the mean of 100 such realizations and the corresponding standard deviation map are shown in Fig. 6. Comparing Figs. 5 and 6, one could observe the displacements of the peak regions in both the rainfall and the standard deviation maps. Compared to Fig. 5, a visible reduction in the standard deviation is observed in Fig. 6, showing the reduced estimation uncertainty by integrating the uncertainty of the radar-indicated rainfall pattern.

Figure 6 Stage 2 results: mean and standard deviation of 100 realizations (mm), obtained using the expected quantile map as the reference.

Figure 7 Stage 3 results: mean and standard deviation of 100 expected realizations (mm).

Stage 3 involves a simulation strategy that is slightly more complicated than before. The simulation is applied independently using the shifted quantile map associated with a positive probability as the reference, and the single simulation cycle is applied to all the components possessing a positive probability. Finally, the expectation of these realizations is computed, termed the expected realization hereafter. Yet the computational cost to obtain such an expected realization is much higher as multiple realizations are required to make up an expected realization (see Table 1 for details). Similarly, the mean of 100 expected realizations and the corresponding standard deviation map are shown in Fig. 7. Comparing Figs. 6 and 7, the positions of the peak regions remain unchanged in both the rainfall and the standard deviation maps. Yet compared to Fig. 6, the reduction in the standard deviation is remarkable in Fig. 7, showing further reduction in terms of the estimation uncertainty.

Figure 8 Stage 3 results: four randomly selected expected realizations (mm).

In addition, the mean field and the standard deviation map of the results from Stage 3 are much smoother compared to the other alternatives shown previously. Figure 8 displays four randomly selected expected realizations to show the similarity of different expected realizations and the smoothness of individual expected realizations.

5.2 Simulation Session II

In the previous session, the objective function only contains the component field pattern. In this session, the component directional asymmetry (abbreviated as asymmetry hereafter) comes into play in the objective function. It should be noted that the evaluation cost of the objective function in Session II is higher than that in Session I, and more time is needed to generate a single realization if the same stopping criterion is used (see Table 1 for details).

We still adopt the three-stage evolution when applying the PA algorithm as in Simulation Session I. Differences between realizations from Session I and Session II do exist, but they are not that remarkable. Therefore, in the presentation of the results from Session II, we omit the results from Stages 1 and 2 and only display the results from Stage 3 in Figs. 9 and 10 (both panels on the right of each figure). The corresponding results from Session I are presented on the left just for the sake of comparison and identifying the effect of adding the component asymmetry in the objective function. Comparing the two mean fields in Fig. 9, the distinction is in the peak values: the two have similar estimates for the peak values of the large rain cell yet different estimates for those of the small cell. As for the comparison of the standard deviation maps in Fig. 10, a slight reduction in the standard deviation of Session II can be observed.

Figure 9 Mean of 100 expected realizations (mm) from Session I (a) and Session II (b).

Figure 10 Standard deviation of 100 expected realizations (mm) from Session I (a) and Session II (b).

Figure 11 (a) Reference asymmetry function. (b, c) Mean and standard deviation of 100 simulated asymmetry functions. Unit: –.

From the results shown in Figs. 9 and 10, the effect of adding the component asymmetry in the objective function seems to be minor, which is due to the fact that the two components (field pattern and asymmetry) share a special relationship: high similarity in terms of the pattern between the reference and the simulated field suggests high similarity in the asymmetry of the two as well. Yet, it does not work in reverse.

To show the capability of the proposed method in terms of fulfilling the component asymmetry, the mean of 100 simulated asymmetry functions is displayed in the middle in Fig. 11 in comparison with the reference asymmetry function (the one evaluated from the reference field) on the left. As for the analysis of the asymmetry function: first, it is a symmetrical function with respect to the origin (0, 0) (see the formula given in Eq. (8)); second, the value of each pixel is obtained by (a) shifting the simulated field by the coordinate of the pixel and then (b) evaluating Eq. (8) for the intersection of the shifted and the original field. This statistic could be computed more efficiently in Fourier space by generalizing the approach presented in Marcotte (1996). As shown in Fig. 11, the extremes of the asymmetry functions for both the reference and the mean simulated ones approach ±0.3. If one computes the asymmetry function for a spatial field obeying a multivariate Gaussian distribution, where symmetric behavior of the field is assumed, the function value approaches zero.

As shown in Fig. 11, the reference and the mean asymmetry are barely distinguishable. The standard deviation map, displayed on the right, reveals the small variability between the simulated asymmetry functions. Note that the presented asymmetry functions are evaluated from realizations of Stage 1. But a similar standard in terms of the fulfillment of the component asymmetry can be achieved by all the other results from Session II.

It is worth mentioning that a trick is used to reduce the computational cost substantially at a fairly low cost in terms of the estimation quality. In Sect. 3, when using multiple shifted quantile fields to produce the expected realization, one should be aware that these individual fields make fairly different contributions to the final estimate. By analyzing the Lorenz curve of the contribution of these individual fields (as shown in Fig. 12), i.e., the weights indicated by the probability matrix, one could see that the lowest several weights contribute very little to the accumulated total. Specifically, the top 9 weights (out of 23) contribute to 90 % of the accumulated total, the top 13 contribute to 95 %, the top 18 contribute to 99 %, and the top 19 contribute to 99.5 %. Being slightly conservative, we have chosen to use the top 19 contributors to produce the expected realization. Note that the top 19 weights should be scaled such that the sum of them equals 1. We have tested this trick (using 19 to represent 23) on the results from Session I (Stage 3). As expected, the difference between the expected realizations, obtained with and without using the trick, is tiny: the maximum difference is 0.050 mm, the minimum difference −0.058 mm, and the mean difference 0.001 mm. The results, shown in Figs. 9 and 10, are produced by using the trick, which helps in saving the computational cost (see Table 1 for details).

Figure 12 Lorenz curve of the contribution of the individual shifted field, where the x axis denotes the cumulative share of population ordered by contribution from the lowest to the highest, and the y axis denotes the cumulative share of contribution.

Table 1 Mean time consumption to generate a realization and an expected realization (with and without using the trick).

The above is based on the performance of a normal laptop.

6 Conclusions

The focus of this paper is to simulate rainfall fields conditioned on the local constraints imposed by the point-wise rain gauge observations and the global constraints imposed by the field measurements from weather radar. The innovation of this work comes in three aspects. The first aspect is the separation of the global and local constraints. The global characteristic of PA makes it a powerful methods for handling the global constraints. Thus PA is only used to deal with the global constraints, and the local ones are handed over to residual kriging. The separation of different constraints makes the best use of PA and avoids its insufficiency in terms of the fulfillment of the local constraints. The second is the extension of the PA algorithm. Except for annealing the system temperature, the number of perturbed phases is also annealed during the simulation process, making the algorithm work more globally in initial phases. The global influence of the perturbation decreases gradually at iterations as the number of perturbed phases decreases. The third is the integration of the uncertainty of the radar-indicated rainfall pattern by (a) simulating using the expectation of multiple shifted fields as the reference or (b) applying the simulation independently using multiple shifted fields as the reference and combining the individual realizations as the final estimate.

The proposed method is used to simulate the rainfall field for a 30 min event. The algorithm of PA is applied using different scenarios: with and without integrating the uncertainty of the radar-indicated rainfall pattern and with different objective functions. The capability of the proposed method in fulfilling the global constraints, both the field pattern and the directional asymmetry, is demonstrated by all the results. Practically, the estimates, obtained by integrating the uncertainty of the radar-indicated rainfall pattern, show a reduced estimation variability. And obvious displacements of the peak regions are observed compared to the estimates, obtained without integrating the uncertainty of the radar-indicated rainfall pattern. As for the two options to integrate the uncertainty of the radar-indicated rainfall pattern, (b) seems to be superior to (a) in terms of the substantial reduction in the estimation variability and the smoothness of the final estimates. As for the two simulation sessions using different objective functions, the impact of adding the component directional asymmetry in the objective function does exist, but it is not that prominent. This is due to the special relationship between the two global constraints: high similarity in the field pattern is a sufficient condition for high similarity in the directional asymmetry function (although the inverse is not true). Yet, compared to the results using the objective function containing solely the component field pattern, a slight reduction in the estimation variability is observed from the results using the objective function containing also the component directional asymmetry.