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 tempo-spatial variability and the intermittent nature of precipitation, it is challenging to obtain QPE with adequate accuracy. This paper aims at simulating rainfall fields 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 the modified phase 5 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 10 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. Copyright statement. The authors retain the copyright of this article.

The impact of directional asymmetry does not necessarily increase with accumulation time (see Figure 1, right). The impact is more dependent on the event, and even instantaneous map of rain rate can present significant asymmetry behavior.
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 30 of geostatistical methods aiming at simulating conditional Gaussian fields with a given covariance function, such as turning bands simulation, LU decomposition-based methods, sequential Gaussian simulation, etc. (see Deutsch et al., 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 non-linearly to the simulated targets, third-or higher-35 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, meta-heuristics method, simulated annealing (SA) (Kirkpatrick et al., 1983;Geman and Geman, 1984;Deutsch, 1992;Deutsch et al., 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 that the perturbation, or 40 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 45

SA.
A remarkable feature of PA is that it is a global method: any perturbation imposed on the phase component is reflected on 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 imparts it a valuable methodology for global constraints. However, PA is found to be insufficient in the treatment of local constraints (Hörning and 50 Bárdossy, 2018). Note that by local constraints, we refer primarily to point equality constraints, when the total number of the 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 in handling the local constraints can be utilized in turn. However, in the general case, the 55 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 work more globally at initial phases of the simulation. The global impact of the perturbation is 60 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 Section 2, including three stages: pre-simulation, simulation, and post-simulation. 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 Section 5, the proposed algorithm is used to simulate the rainfall field for a 30-min-event, where 65 different simulation scenarios are applied. Section 6 ends this paper with conclusions. Figure 1 summarizes the procedure of simulating surface rainfall fields using the algorithm of PA, including three stages:

Methodology
pre-simulation (PreSim), simulation (Sim) and post-simulation (PostSim). Each stage and the corresponding sub-stages are described in the following subsections in the same sequence as shown in the flowchart.  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.

75
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 1  ···  r K and u 1  ···  u K .

80
(b) Set the spatial intermittency u 0 , i.e. the quantile corresponding to zero precipitation. u 0 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 1 to estimate u 0 : where r 1 is the smallest gauge observation corresponding to u 1 . As an alternative, if the small quantiles of the field are 85 not properly sampled by the rain-gauges, we propose to estimate u 0 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 : with r k 1 , r k being the two nearest neighbors of r (r k 1  r  r k ), and u k 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 Equation 3. We have learned from practice that the exponential extrapolation tends to 95 over-estimate the rainfall extremes, so a linear component is used to restrict the extrapolation result.
where, the only parameter of the exponential distribution is determined from the last pair (r K , u K ): It is recommended that the methodology described above is used to estimate the distribution function of accumulated rainfall, 100 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 105 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 Section 4, is necessary to maintain the accuracy of radar-indicated spatial ranks.

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 where 1 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 Section 2.1.1.

115
We impose two kinds of constraints: local and global. As has been explained in the introduction, PA is a powerful method to handle 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 120 simulated Gaussian field and are compared with the radar-based Gaussian field Z ? , as obtained in Equation 6. And we term Z ? the reference field, hereafter.
The first global constraint, field pattern, requires that the simulated field should be similar to the reference field. The similarity of the two fields is quantified by the Pearson correlation coefficient: In the ideal case, ⇢ Z,Z ? equals 1, and we use the difference, (1 ⇢ Z,Z ? ), to measure the distance from the ideal.

specified it is a CDF not a PDF
The second global constraint is directional asymmetry, as given by 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 ⇡ h and is the cumulative standard normal distribution function. Directional 130 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 1 norm, kA Z A Z ? k L 1 .

135
Different schemes could be used to combine the two components, linear or non-linear, 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: where w is the relative weight of the component directional asymmetry and A scal , as expressed in Equation 10, is the scaling 140 factor that scales the L 1 norm of the difference of the two asymmetry functions between 0 and 1.
It is worth mentioning that the proposed methodology is flexible and any global statistic could be used to constrain the simulated fields, and one could even combine several statistics in the objective function.

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.

150
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 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 non-monotonic adaptive cooling, have been made by, e.g., Nourani and Andresen (1998), Martín and Sierra (2009) and others. The results show that annealing works properly when the cooling curve has a moderate slope at the 160 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 0 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 165 corresponding acceptance rate and improvement rate are computed. The acceptance means the perturbation is being accepted by the system and the system state is being updated, whether the perturbation brings improvement to the system; while the improvement means the perturbation does decrease the objective function value. In short, a perturbation bringing improvement must be an accepted perturbation, yet an accepted perturbation is not necessarily bringing improvement to the system. We refer to the experiment as a temperature cycle, where N perturbations are being implemented and the two mentioned rates are being 170 computed at a fixed temperature.
T 0 is first set by starting a temperature cycle with an initial guess of T 0 . If the acceptance rate is lower than the predefined limit, say 98%, then increase the temperature, and vice versa. The goal is to find a T 0 with a relatively high acceptance rate.
Then, T min is set by decreasing T 0 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 175 it is again a compromise between a satisfactory destination and the computational cost.
We note that from the determination of T 0 , the total number of temperature cycles m (to anneal T 0 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 0 , T min and L. The temperature at Iteration l is computed as where ↵ T is the attenuation factor of the system temperature, given as 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 0 (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 where ↵ N is the attenuation factor with respect to the number of perturbed Fourier phases. It is computed as The logic behind annealing the number perturbed phases is to render the perturbation have a more global impact initially when the distance from the destination is relatively large, and decrease the impact of the perturbation gradually as the opti-190 mization process develops.

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 band simulation (Journel, 1974), or the Cholesky transformation of the covariance matrix.

195
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. zero. If covariance models with asymptotic ranges (e.g. exponential, Matérn, Gaussian covariances, etc) are employed, the extension in domain size could be significant (Chilès and Delfiner, 2012).

PA Main Cycle
The starting point of PA is a Gaussian random field Z K with the prescribed spatial covariance, as described in Section 2.2.2.
Here a little modification, in our case, is that the values at the observational locations are fixed by residual kriging, as indicated 205 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  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 0 ] and [1, · · · , V 0 ], respectively, where U 0 , V 0 are the highest frequencies considered for perturbation in x-and y-direction, respectively. Note that zero frequencies in both directions should be excluded from the perturbation.
One could use the entire frequency range or select a sub-area to impose the perturbation. 4. N l phases ✓ n are randomly drawn from the uniform distribution, [ ⇡, ⇡). 215 5. The Fourier coefficients at the selected locations (u n , v n ) and the corresponding symmetrical locations are expressed as Z K [u n , v n ] = a n + jb n (16) where U, V are the number of grid points in x-and y-direction, respectively (also the highest frequencies in both directions). These coefficients are then updated in terms of the Fourier phases using ✓ n ; while the corresponding amplitudes 220 remain unchanged: The perturbed Fourier transform is denoted asZ. The corresponding perturbed spatial random field is obtained by the inverse DFT: 6. Due to the perturbation,Z no longer satisfies the point equality constraints (note the removal of the subscript "K"). Thus kriging is applied for the residuals z k Z (x k ), where z k are the point equality constraints defined in Equation 5. The results of kriging, r, are superimposed onZ: As kriging is a geostatistical method that depends only on the configuration of the data points, the weight matrix of individual grid point 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.Z K is then subject to the objective function defined in Equation 9. If O(Z K ) < O(Z K ), the perturbation is accepted.
Otherwise, the perturbation is accepted with the probability: 8. If the perturbation is accepted, the system state is updated, namely Z K , Z K and O(Z K ) are replaced byZ K ,Z K and O(Z K ), respectively.
9. If the stopping criterion is met, stop the optimization and Z K is a realization in Gaussian space, satisfying the predefined optimization criterion; otherwise, go to Step 2.

Post-simulation
The Gaussian field Z K is transformed into a rainfall field by: where is the cumulative standard normal distribution function and G 1 is the inverse of the distribution function of surface rainfall, as obtained in Section 2.1.1. The resultant R is a realization of the surface rainfall field. 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 to 1, the more accordance (or less conflict) of radar and gauge data.

250
The conflict of the two could be partially explained by the fact that weather radar is measuring 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 255 and Bárdossy (2019) is adopted. Here we recall the procedure briefly in three steps.
(1) A rank correlation matrix ⇢ S ij 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 .

260
Both h ij and ⇢ S ij have the same spatial resolution as U . The center of h ij is (0, 0), with the corresponding entry in ⇢ S ij denoted as ⇢ S 00 , meaning that zero shift is imposed on U . The farther a grid cell is from the center, the larger shift is imposed on U . 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 ⇢ S ij : where H is a monotonic function, such as power function, exponential function or logarithmic function, etc. The first derivative (dH/dx) matters, as a large dH/dx means that more weights are assigned to those displacement vectors (or shifted fields) producing a higher rank correlation. Our choice of function H is simply: H(x) = x 2 .
(3) The probability matrix P ij is scaled to ensure that the sum of all entries equals 1: 270 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 increasing 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 19km ⇥ 19km (see Section 4 for details), the selection of the extent is limited by the gauge data 275 availability. With this spatial extent, quite often we could 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 40km ⇥ 40km), 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.

280
The first option: the expectation of the shifted quantile fields is computed using the probability matrix, and then the corresponding Gaussian field is computed, as given in Equation 26. The result is used as the reference field when applying the PA algorithm.
The second option: those (marginal-converted) shifted quantile fields, 1 (U ij ), bearing a positive probability are consid-285 ered as reference fields and the algorithm described in Section 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: 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 Section 5: Application.

Study Domain and Data
The study domain is located in Baden-Würtemberg in the southwest of Germany. As shown in Figure 2, it is a square domain with the side length of 19km. The domain is discretized to a 39 ⇥ 39 grid with the spatial resolution of 500m ⇥ 500m. A gauge 295 network consisting of 12 pluviometers is available within the domain, as denoted by the red dots in Figure 2.
The radar data used in this study is the raw data with 5min temporal resolution, measured by Radar Türkheim, a C-band radar operated by the German Weather Service (DWD). Radar Türkheim is located about 45km from the domain center, as denoted by the red star in Figure 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) re-projection 300 from polar to Cartesian coordinates and (4) clip the square data for the study domain. All these quality control steps are operated under the environment wradlib, an open-source library for weather radar data processing (Heistermann et al., 2013).

Application
We select a 30-min event to apply the algorithm of PA. The event is selected not only due to the relatively prominent rain intensity reflected by the rain-gauge data, but more importantly, the event is properly captured by a few rain-gauges unevenly 305 distributed in the domain of interest, as shown by the red dots on the left panel in Figure 3. From the figure, it is 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 are more or less evenly distributed. One might have noticed that the smallest sampled quantile is 0.53. Yet, in this case, u 0 = 0.26 (the spatial intermittency) and the quantile value 0.53 corresponds to the rainfall value 0.28mm after the re-ordering.

310
The conflict of radar and gauge data is obviously reflected in Figure  However, using the algorithm described in Section 3 to evaluate the uncertainty of the radar-indicated rainfall pattern, one 315 can obtain multiple distribution functions of surface rainfall, as shown by the grey lines on the right panel in Figure 3. 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 on the right panel in Figure 4. Also with 320 the probability matrix, one can compute the expectation of the shifted quantile fields, as shown on the left panel in Figure 4.
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 Figure 3 and 4.
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 in the middle. 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 to 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 guaranteed.
We open up two simulation sessions, depending on the different objective functions used when applying the PA algorithm.

330
In the first session, the objective function contains only the component field pattern: O(Z) = 1 ⇢ Z,Z ? . In the second session, the component directional asymmetry comes into play and the objective function expressed in Equation 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.

Simulation Session I 335
We present an evolution in terms of the simulation strategy, where the algorithm of PA is applied differently in 3 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 Section 3; (3) simulating independently using those shifted quantile maps with a positive probability as the reference and computing 340 the expectation, i.e., integrating the uncertainty of the radar-indicated rainfall pattern via Option 2 in Section 3.
Stage 1, simulating rainfall fields, using the original quantile map as the reference, means the uncertainty of the radarindicated 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 Figure 5 (right), which shows the fulfillment 345 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 radarindicated rainfall pattern via the first option, as described in Section 3. Similarly, the mean of 100 such realizations and the corresponding standard deviation map are shown in Figure 6. Comparing Figure 5 and 6, one could observe the displacements 350 of the peak-regions in both the rainfall and the standard deviation maps. Compared to Figure 5, a visible reduction in the standard deviation is observed in Figure 6, showing the reduced estimation uncertainty by integrating the uncertainty of the radar-indicated rainfall pattern. 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 355 to all the components possessing a positive probability. Finally, the expectation of these realizations is computed and 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 Figure 7. Comparing Figure 6 and 7, the positions of the peak-regions remain unchanged in both the rainfall and the standard deviation maps. Yet compared to Figure 6, the reduction 360 in the standard deviation is remarkable in Figure 7, showing the further reduction in terms of the estimation uncertainty.
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

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 so does the time 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. The differences between 370 realizations from Session I and Session II do exist, but are not that remarkable. Therefore, in the presentation of the results from Session II, we omit the results from Stage 1 and 2 and only display the results from Stage 3 in Figure 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 Figure 9, the distinction consists in the peak-values: the two have similar estimates for the peak-values of the large rain cell, 375 yet different estimates for that of the small cell. As for the comparison of the standard deviation maps in Figure 10, a slight reduction in the standard deviation of Session II can be observed. From the results shown in Figure 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 380 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 Figure 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 Equation 8; second, the value of each pixel is obtained by (a) shifting  in Marcotte (1996). As shown in Figure 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.

390
As shown in Figure 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 played to reduce the computational cost substantially at a fairly low cost of the estimation 395 quality. In Section 3, when using multiple shifted quantile fields to produce the expected realization, one should be aware that these individuals have fairly different contribution to the final estimate. By analyzing the Lorenz curve of the contribution of these individuals (as shown in Figure 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 400 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 playing the trick, is tiny: with the maximum difference 0.050mm; the minimum difference 0.058mm and the mean difference 0.001mm.
The results, shown in Figure 9 and 10, are produced by using the trick and it helps in saving the computational cost, see Table 1 405 for details.  The above is based on the performance of a normal laptop.

Conclusions
The focus of this paper is to simulate rainfall fields conditioned on the local constraints imposed by the point-wise raingauge observations and the global constraints imposed by the field measurements from weather radar. The innovation of this work comes in three aspects. First, separate the global and local constraints. The global characteristic of PA imparts it a 410 powerful methods in 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. Second, extend 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 at initial phases. The global influence of the perturbation decreases gradually at iterations as the number 415 of perturbed phases decreases. Third, integrate 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; with different objective functions.

420
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 425 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 is not that prominent. This is due to the special relationship between the two global constraints: high similarity in the field pattern is the 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 430 the estimation variability is observed from the results using the objective function containing also the component directional asymmetry.