NowPrecip: localized precipitation nowcasting in the complex terrain of Switzerland

The new precipitation nowcasting system “NowPrecip” is introduced and its methodology is described in detail. NowPrecip is an area‐tracking application, able to generate ensembles of precipitation nowcasts. It supports the paradigm of seamless forecasting: nowcasts start from the most recent radar observation and evolve smoothly, through merging, towards the numerical weather prediction (NWP) ensemble members. NowPrecip relies on the following three pillars. (a) A geostatistics‐based optical flow algorithm, named “NowTrack”, able to determine the scales of interest automatically and equipped with an outlier control scheme. Such a scheme is particularly useful in capturing situations like fine‐scale rotation accurately. (b) A localized architecture motivated by the need to address scenarios where distinct rainfall patterns characterize different parts of the domain. This is mostly relevant for complex terrains like Switzerland, but is also useful for areal nowcasting operating on extensive regions. (c) Techniques to compute and incorporate localized growth and decay into the computation of nowcasts. This is especially useful for Alpine terrain, where there is often systematic growth and decay associated with orographic features. The growth and decay factors are based on analysis of the available NWP product, allowing for a more efficient coupling between the nowcasting sequence and the numerical model output. The main scope of this work is to present the methodology of NowPrecip in detail, but a verification of several representative events is also provided, using both deterministic and probabilistic measures.

detects and isolates precipitating cells, treats and tracks them as individual entities, and forecasts their evolution deterministically or probabilistically. It focuses on severe storms and produces forecasts for lead times in a typical range between 15 and 60 min. Area tracking, on the other hand, aims to track and extrapolate a continuous precipitation field that expands over the entire available domain (Hilst and Russo, 1960;Wilson, 1966;Bellon and Austin, 1978;Golding, 1998;Germann and Zawadzki, 2002;Bowler et al., 2006;Haiden et al., 2011). It provides predictions for several hours (6 hr is typical) and nowadays gradually merges the nowcast with the corresponding numerical weather prediction (NWP) forecast. The skill of such merged nowcasts is generally superior to those of the NWP alone or the nowcast alone for at least a few hours (Lin et al., 2005;Mandapaka et al., 2012). Area-tracking algorithms consist of two main parts. First, the observed precipitation is analyzed and a motion field is generated to describe the apparent radar echo motion. Second, this motion field is used to extrapolate the most recent rainfall observation to short lead times. In advanced schemes, the rainfall evolution of the forecast field at different scales is assessed using autoregressive schemes that include spatially correlated noise in the form of a shock term. This term is consistent with the spatial statistical characteristics of the most recent observations and ensures that the lifetimes of features at different scales are reproduced correctly (Seed, 2003;Bowler et al., 2006).
It is important to forecast precipitation accurately for short lead times, since rainfall has a significant social and economic impact. Heavy rainfalls often function as trigger events for geographically localized hazards like flash floods, debris flows, and landslides (Sättele et al., 2012). Such hazards are particularly prominent in the complex Swiss terrain, where the Alps cover 60% of the country and the landscape is covered by numerous rivers and over 1,500 lakes (Spreafico and Weingartner, 2005). Moreover, severe rainfall is often associated with a high density of lightning flashes, which are responsible for causing damage to power supplies. Precipitation is also connected to strong wind gusts and shear turbulence, which affect aviation severely. Extreme rainfall causes, directly or indirectly, human injuries or fatalities, as well as damage to material assets, traffic lines, infrastructure, forestry, and agricultural land (Hilker et al., 2009). One hundred and two fatalities were reported in Switzerland for the period 1972-2007, while the cost of damage caused by thunderstorms alone was estimated to be 1.85 billion Euros. Thunderstorm-related damage comprises about 23% of the total cost of storm-related damage for this period (75% or 6 billion Euros has been attributed to long-lasting rainfall). 95% of the total cost was incurred between May and October. Warnings, based directly on nowcasting predictions or using nowcasting outputs as intermediate information, for example, forcing for hydrological models, are critical in reducing both human fatalities and financial costs. This is because, in most cases, the impact is not only "weather-sensitive" but also "weather-information-sensitive", meaning that correct information can minimize the impact (Cuo et al., 2011). Even when precipitation is not intense, accurate nowcasting has a positive financial and social impact in the sectors of transportation, tourism, business, and daily activities.
Nowcasting relies on the concept of "persistence": precipitation patterns stay relatively unaltered for long enough to be predicted if extrapolated into the future. A radar two-dimensional image is essentially a horizontal projection of an underlying vertical process, a continuous and turbulent physical flow of mass and energy among a wide range of scales. At large scales, patterns maintain their coherence longer, while small features are short-lived. This was an observation that, in part, defined the zeitgeist in the early history of nowcasting and naturally raised expectations for predicting the future evolution of radar echoes. There is a plethora of relevant literature from the 1960s (Noel and Fleisher, 1960;Boucher and Wexler, 1961;Wilson and Kessler III, 1963;Wilson, 1966), and probably even earlier, in works that are difficult to access today (Blackmer, 1953). The thrust of this early research led to the work of Zawadzki (1973), who refined the concept considerably by studying the variability of rainfall rate as a function of space and time separately, and has stayed active in this area through numerous contributions since (Germann and Zawadzki, 2002;Turner et al., 2004;Germann et al., 2006a;Radhakrishna et al., 2012).
Today extrapolation-based forecasts are still superior to NWP forecasts, at least for short lead times (0-3 hr), and there are several reasons for that (Simonin et al., 2017). First, NWP schemes employ assimilation processes (the merging of observations distributed in time with the dynamical numerical model of the flow) in order to determine the state of the atmosphere as accurately as possible and produce appropriate initial conditions (Talagrand, 1997). Even for high-resolution models, these processes have limitations that can result in discontinuities between their outcome and the most recent radar precipitation observations (Reyniers, 2008). Second, NWP forecasts take time to be completed: they are typically ready well after real time, even for rapid update cycles (Pierce et al., 2012). Consequently, small-scale weather features such as thunderstorms, with sizes up to a few tens of kilometres and lifetimes between 30 and 60 min, cannot yet be forecast by NWP with sufficient spatial and temporal accuracy. A third reason concerns the statistical characteristics of precipitation. Rainfall is characterized by strong spatiotemporal variability, due to the complexity of the underlying macro-and microphysical processes on scales smaller than the mesoscale or even comparable to or smaller than the grid scale of NWP models. As a consequence, accurate fine-scale predictions by numerical models are challenging, especially when convective and orographic processes play a role (Golding, 2000;Ebert et al., 2003;Cuo et al., 2011). Problems with displacement errors have been reported to be the main source of quantitative precipitation forecasting (QPF) errors (Ebert and McBride, 2000). Such errors affect hydrological applications, especially streamflow forecasting operations (Cuo et al., 2011). Finally, precipitation is not a prognostic variable in NWP models and can only be assimilated indirectly through other variables, for example latent heat nudging.
Nowadays the expectations for a high-end nowcasting scheme are high. Preferable characteristics of a nowcasting system are as follows. (a) The nowcasting system should be probabilistic, that is, able to generate a multimember ensemble output, where each member is consistent with the statistical properties of the most recent observations. The spread of the ensemble should capture the uncertainty in hand adequately and ideally should include possible extremes, like heavy rainfall cells. (b) It should blend seamlessly with the members of an available NWP product (nowcasting is superior for short lead times, until the initial conditions wash out, then NWP takes over). (c) It is preferable to have a localized architecture, that is, to be able to describe non-stationary scenarios, characterized by diverse behaviours in different parts of the domain (a situation that occurs when the domain size is large or the topography is complex), adequately. (d) It should analyze and employ information about current growth and decay, especially the persistent kind, which is usually associated with topography. (e) It is desirable to rely on a motion field that represents the motion accurately at all scales, from large to small, that play a significant role. Moreover, the field produced has to be continuous, smooth, and artefact-free, even in regions with zero precipitation. (f) A final, relatively complex, expectation is to incorporate information from object-based nowcasting into area-based nowcasting. For example, satellites may occasionally detect convective initiation minutes before its signature is detected by radar. This is important for alerts, and ideally a nowcasting algorithm should be able to incorporate this information by seamlessly introducing a convective cell into the evolution .
This article presents a set of new nowcasting methods in three main areas: (a) optical flow, (b) localization, and (c) growth and decay. These components are integrated into the novel MeteoSwiss nowcasting system called "NowPrecip", a system that has been developed by R Development Core Team (2016) and belongs to the same suite of applications as the raingauge-radar combination system of MeteoSwiss, "CombiPrecip" (Sideris et al., 2014). The algorithms can also be seen as standalone techniques and can be employed independently in other nowcasting systems. In order to make the methodology beneficial and reproducible by others, the authors have opted for presenting it in relative detail. Therefore, the emphasis of this article is dedicated to the demonstration of the algorithms and associated examples. This is accompanied by a verification section, which is shorter in length, but is based on a sufficiently sizeable set of events to serve as a first solid assessment of the new system. Our focus in this article is therefore clearly on the methodology rather than extensive verification of the material presented. Such verification should be presented in future work.
The article is organized as follows. In section 2 the data and the list of precipitation events are presented. In section 3 the methodology is presented. This section is divided into three algorithms: the optical flow, the core localized architecture, and the growth and decay algorithm. In section 4 the verification is presented, followed by our conclusions in section 5, and a description of future work in section 6.

DATA
Five C-band dual polarization radars are currently operated in Switzerland by MeteoSwiss. Their precipitation estimates are subject to detailed quality control and produce a two-dimensional precipitation composite in a domain that covers Switzerland and parts of the neighbouring countries. The radar domain extends over a 710 × 640 km Cartesian grid, and has spatial and temporal resolution of 1 km and 5 min respectively. The radar composite is refined further by combining it with raingauge measurements to produce a merged product called Com-biPrecip (Sideris et al., 2014), but for the needs of this article the radar-alone composite has been used. Information on the Swiss weather radar network and the steps for quantitative precipitation estimation and quality control have been reported in the literature (Germann et al., 2006b;. The merging process of NowPrecip uses the outcome of the NWP system COSMO-CH (Schättler et al., 2008;Baldauf et al., 2011). In particular, COSMO products from the versions COSMO-1 (deterministic, resolution 1 km) and COSMO-E (ensemble of 21 members, spatial resolution 2.2 km) are employed. For the purposes of NowPrecip, the COSMO products were cropped to span the aforementioned radar domain, while COSMO-E was downscaled to a spatial resolution of 1 km (the same as the radar product). The temporal resolution of COSMO products is 10 min. The output of NowPrecip also has a temporal resolution of 10 min and spatial resolution of 1 km. In the operational environment, COSMO-1 has an update cycle of 3 hr while COSMO-E has an update cycle of 12 hr. Sixteen events were chosen manually to demonstrate the methodology and capabilities of NowPrecip. Included are both stratiform events, which happen mainly during winter, and pure convective events characterized by isolated distinct cells, which happen mainly during the summer months. Several events are characterized by convective precipitation embedded in a stratiform background, a common and occasionally dangerous type of precipitation in Switzerland during summer. In some of the events, orographic effects are dominant; in others, small-scale rotation takes place. These events comprise a representative sample of types that are common in Switzerland and several of them can be characterized as rather challenging for a nowcasting system. A list of the events and their main characteristics can be seen in Table 1.

Optical flow
The sequence of the two most recent observed rainfall fields can be represented by a motion field. A geostatistics-based algorithm, coined as "NowTrack", is introduced here to estimate this motion field. This algorithm is comprised of three parts. (a) Tracking of precipitation using a hierarchical cross-correlation technique for pattern matching of a small number of pixels (Rinehart and Garvey, 1978;Anandan, 1989). These velocity estimates consist of the "measurements". (b) A quality-control process to exclude measurements that are inconsistent with their neighbours. (c) Use of these remaining measurements to perform kriging to estimate the velocity fields at every pixel of the rainfall map. These steps are now described in detail.
1. Two successive radar rainfall observations, at times t 1 and t 2 = t 1 + , are used for the estimation of the velocity field. The observation at time t 2 is the most recent available radar observation. The two observations represent 10-min precipitation accumulations in mm units, and = 10 min. N p = 300 locations are selected randomly for the precipitation observation at time t 2 . These locations have coordinates ( The precipitation at each of these locations must be larger than zero. A box centred on each location (x i , y i ) and extending 48 pixels at each direction (i.e., the box size is 97 × 97 pixels) is defined. This choice of box size was found to be the best after examination of errors in several examples. 2. This box is matched (in a maximum cross-correlation sense) with precipitation in a corresponding box in the time frame t 1 . In the observation at time t 1 , a search region is defined, centred at the location (x i , y i ) and with an extent of 18 pixels in each direction. The choice of 18 pixels is adequate for typical meteorological situations in Switzerland, since the time difference between two consecutive observations is = 10 min = 1∕6 hr, which implies that the maximum speed of the precipitation patterns that can be captured with the proposed scheme is 6 × 18 = 108 km⋅hr −1 (which did not occur in this study). Now the 97 × 97 pixel box centred at (x i , y i ) at time frame t 2 is compared with the moving 97 × 97 pixel box centred at (x i ± s 0 , y i ± s 0 ) at time frame t 1 , where s 0 ∈ {0, 6, 12, 18}. In particular, the Pearson correlation coefficient between the original box at t 2 and the moving box at t 1 is computed. This search takes place every 6 pixels for computational efficiency, because the goal in this first search is to find the approximate location where the best match takes place. 3. Through this process, the location of the highest correlation in the 6-pixel resolution grid is established. Let this location be (x j , y j ). This process is now repeated for a box with size 49 × 49 pixels (half the extent of the one in step 2), but now the search takes place in the region centred at (x j ± s 1 , y j ± s 1 ) with s 1 ∈ {0, 3, 6, 9} and within a box region that extends 9 pixels in each direction. In this hierarchical way, the coordinates of the maximum correlation between the two frames are established more carefully. Let the new best-match location have coordinates (x k , y k ). 4. The process is repeated for a third and final time, now for a box that has size 25 × 25 pixels (half extent of the one in step 3), but now the search takes place in the region centred at (x k ± s 2 , y k ± s 2 ) with s 2 ∈ {0, 1, 2, 3} and within a small box region that extends 3 pixels in each direction. Let the final best-match location have coordinates (x l , y l ). This hierarchical pattern-matching process is depicted for two examples in Figures 1 and 2. In the large panel on the left, one can see the search boxes, one inside the other, depicted in white dashed lines, while the black squares represent the matching-pattern boxes. The two small panels on the right show the two 97 × 97 pixel boxes at times t 1 and t 2 centred on the same Lagrangian frame (the frame that moves with the storm) for comparison. This centring is based on the motion vector computed as mentioned above. In practice, if the derived motion vector were incorrect, then an obvious spatial misplacement would be evident between the two panels.
This algorithm returns the best-match location, therefore the two components of the displacement, Δx and Δy, for each measurement, and the corresponding Pearson correlation coefficient between the final 25 × 25 pixel boxes. It is possible to opt for different box sizes in this scheme. However, the algorithm as described here was able to capture the large majority of the features of importance for the needs of areal tracking. Similar hierarchical approaches have been presented before (Anandan, 1989). Hierarchical methods are less computationally expensive and minimize the chance for accidental mismatches by guiding the matching process, from coarse matches to finer ones (Beauchemin and Barron, 1995). 5. The (N p = 300) measurements go through a quality-control process. Some of the measurements are erroneous. Growth and decay occasionally mislead the tracking processes and cause errors. An outlier control scheme is applied to exclude the one third of the measurements that have the largest discrepancies compared with their neighbours. This is a simple but robust condition, motivated by the operational character of NowPrecip. This process "cleans" the data and results in 200 trustworthy measurements, which will be used later in the kriging interpolation. Consistency checks have been used before, for example, as a component of the nowcasting system short-range warning of intense rainstorms in localized systems (Li et al., 2000), although no details were provided in that particular report. The outlier control operates in two dimensions x and y and is computationally fast. The main question asked is to what extent each measurement is consistent with its neighbouring measurements. If there is a significant discrepancy, then this measurement is excluded. A straightforward solution to this problem is to locate a number of measurements within the vicinity of each measurement, fit these data in a two-dimensional distribution, and compute how many standard deviations away from the central value the measurement in question is. When the number of standard deviations is larger than a certain threshold, the measurement is rejected. This is exactly the idea behind the Mahalanobis distance (Mahalanobis, 1936), which is defined as where Δr = (Δx, Δy) T is the measurement in question, μ = ( Δx , Δy ) T the vector of mean values of the neighbourhood measurements, and S their covariance matrix. The Mahalanobis distance is typically used for outlier detection and essentially checks whether an observation belongs to a distribution, by providing a measure of the observation's distance from the mean value. The covariance function of the neighbourhood measurements determines an ellipsoid, which represents the probability distribution of the set (multinormality The location of the origin becomes specified progressively better in three iterations. The black boxes represent the regions that have to be best-matched between two consecutive time frames at times t and t + . The search begins in the large white box, where the large black box is translated, attempting to locate the position of maximum cross-correlation between the two consecutive time frames. The process repeats twice for smaller search and matching boxes. The last search provides the optimal coordinates of origin (start of the red arrow) for the destination point (end of red arrow) in question. Right panels: boxes of size 97 × 97 pixels centred on the same Lagrangian frame for comparison. The top right panel has been translated using the motion vector computed by the aforementioned hierarchical search system. Any error in the matching process would result in the top right panel image appearing shifted compared with the bottom right panel. The central cross in the two panels has been added to facilitate comparison is assumed). This neighbourhood is a circumference around the location of the measurement which starts at radius 10 km and increases if needed, until it includes at least N p /10 neighbours or the radius reaches a maximum allowed radius of 60 km. The Mahalanobis distance is essentially the distance of the test point from the centre of gravity divided by the thickness of the distribution at that location. After the first pass of this quality control, those measurements found to be more than a threshold of three standard deviations away from the neighbour's set mean value are excluded, and the process is repeated iteratively by decreasing this threshold slightly until 100 measurements are excluded. The remaining 200 measurements are considered to be of high quality and are used as input of a kriging interpolation over the Swiss radar domain.
The original 300 motion vectors computed using the hierarchical algorithm described can be seen in the left panel of Figure 3. It can be seen that some of the arrows inside the three circles are inconsistent with their neighbours. The outcome of the quality-control process can be seen in the right panel of Figure 3. Measurements inconsistent with their neighbours have now been rejected automatically. Figure 4 shows two examples of Mahalanobis outlier control. In the left panel, the measurement (blue triangle) is consistent with its neighbours and obviously belongs to their distribution. Its Mahalanobis distance is small, so this measurement is not excluded. In the right panel, there is a measurement that does not fit well in the distribution of its neighbours. Its Mahalanobis distance is larger than the prescribed threshold and this measurement is excluded. 6. The next step is to produce two fields, for each of the two components of the displacement Δx and Δy on a grid with spacing 1 km. For computational efficiency, this process takes place in two steps. The kriging

F I G U R E 4
Two examples of quality control. In the left panel, the Δx and Δy of the test point (blue) are similar to those in its neighbourhood. The point is obviously located well within the ellipsoid that represents the local distribution. Its Mahalanobis distance is about 0.43, which suggests that this point should not be excluded as an outlier. In the right panel, the test point is far off the main part of the distribution and its Mahalanobis distance, which is about 5.9, reflects that. This point will be excluded as an outlier F I G U R E 5 Motion field semivariograms for Δx and Δy for the example presented in Figure 3 semivariograms for the two components of the displacement are computed and fitted using the exponential model (Goovaerts, 1997). The kriging semivariograms allow for nonzero nuggets to take care of possible microvariability errors in our measurements. Examples are depicted in Figure 5. Two kriging interpolations for Δx and Δy are run independently for a grid with grid length equal to 5 km. This way, the originally irregularly placed measurements are used to generate a regular grid. As a second step, the results of these kriging interpolations go through bilinear interpolations (again one for Δx and one for Δy), which produce the final fields Δx and Δy on the fine grid with grid length equal to 1 km. The final motion fields for six different cases can be seen in Figures 6 and 7. The computational time of this optical flow algorithm increases linearly with the number of measurements. About 15 s are needed for completion of the entire process (pattern matching, outlier control, and interpolation) over the domain of the Swiss radar composite (710 × 640 km). Use of kriging to derive the motion field comes with an advantage, since, apart from estimation of the expected value for each grid point, the solution of the kriging system also provides, at no extra cost, the associated variance ( Figure 8). This provides an idea regarding the estimation error over the domain.
For example, the vectors far away from the nonzero precipitation area are intuitively expected to be quite uncertain, because tracking cannot take place in these regions. The variance provides a quantitative measure of this intuition. This advantage is particularly useful to produce perturbations for the estimated motion field and turn it into probability estimates, that is, equip each member of the ensemble with a slightly different version of the motion field.
Emphasis is now given to two points.
1. A geostatistical interpolation for estimation of velocity fields does not have a need for smoothing constraints similar to those used in other optical flow approaches. This is because kriging is based on fitting a positive-definite model of spatial variability on the experimental semivariogram and relies on this smooth function model to produce the estimation (Goovaerts, 1997). As a result, the estimated fields are characterized by continuity, abrupt changes from pixel to pixel are absent, and artefacts in zero-precipitation regions are minimized. 2. There is no need for a cutoff minimum scale to be prescribed arbitrarily in advance. The smoothness of the estimated motion field relies on the parameter values of the estimated semivariogram, which depends on the input measurements. If the motion field of an event is very smooth, for example, stratiform precipitation, large-scale motion is adequate to describe the motion. On the other hand, cases with local rotation are characterized by much smaller scales and a rougher motion field. Observe that the estimated range of the semivariogram reflects the global decorrelation distance observed within the the entire raster, while local behaviours may vary to different extents. However, the estimated motion field also depends on the measurements themselves. Notice that each measurement is computed independently of each other one, so the scheme can be described as a block estimation of measurements followed by a kriging process, which operates as a smoothing function. This can profit the resulting motion field on certain occasions. Consider, for example, a very turbulent scenario, where two relatively close-by points are characterized by very different motion vectors. In such a scenario, NowTrack should generate two very different vectors for these nearby locations. The kriging semivariogram will be characterized by a small range, suggesting that the variability in this case is large. Therefore, the high velocity field variability will be determined automatically without, for example, the need for a prescribed minimum cutoff scale, through the combination of the fitted semivariogram and block estimation process.
Details of the motion field can be seen in Figure 7. In these images, one sees that NowTrack can handle equally well smooth large-scale motions (as in the top left or middle right panel), cases with orographic effects (as in the top right panel), and locally complex situations like the middle left panel and bottom two panels.

Localized nowcasting architecture
The extrapolation mechanism of NowPrecip relies on an iterative autoregressive scheme. Similar techniques have been presented before (Pegram and Clothier, 2001a;Bowler et al., 2006;Berenguer et al., 2011). NowPrecip builds on these approaches by introducing a localized infrastructure in order to address precipitation evolution in the complex Alpine environment more accurately. Localization is important when there is a significant variation of the statistical properties of precipitation over the domain of interest, for instance, due to orography or because the domain is too large to be characterized by a single precipitation type. Moreover, localization provides an elegant architecture to compute growth and decay and include machine-learning techniques into nowcasting.
The scheme of NowPrecip can be seen in Figure 9. The steps are presented in detail here. The routine for the growth and decay and that for the merging with NWP are presented later in dedicated sections.
1. The six most recent radar observations are retrieved.
As mentioned before, these observations have spatial resolution 1 km and temporal resolution 5 min.
Hereafter, these observations will be referred to as O 0 fields. 2. The six 5-min O 0 radar observations are aggregated into three 10-min accumulations in order to be consistent with the available accumulations of the NWP output. Henceforth, these new 10-min accumulations will be referred to as O 1 fields. The temporal resolution of the O 1 fields also determines the lead-time unit = 10 min of the nowcasting system. We also define T max = 360 min = 6 hr as the maximum lead time chosen for NowPrecip and N = T max ∕ = 36 as the number of frames of the nowcast sequence generated. 3. The motion field is computed using NowTrack, and it is used to advect the fields O 1 (t = − ) and O 1 (t = −2 ) into the frame of the field O 1 (t = 0). Therefore, all three O 1 fields are now "centred" in the same Lagrangian frame. These new "centred" fields will be referred to as O 2 . 4. A normal quantile transformation is applied to turn the nonzero value pixels of the O 2 fields into N(0, 1) distributions, which will be used later in an autoregressive AR(2) scheme (Priestley, 1981;Wei, 2006). In this transformation, only pixels with the nonzero non-NA values are considered. Pixels of the O 2 fields that have zero or NA values are then assigned the minimum value of the transformed distribution, which is typically below −3. This transformation produces a sequence of three frames, which will be referred to as O 3 . 5. Using the O 3 fields, the Yule-Walker coefficients (Yule, 1927;Walker, 1931;Priestley, 1981;Wei, 2006;Clothier and Pegram, 2002;Berenguer et al., 2011) are computed to be used later in the autoregressive scheme. The Yule-Walker coefficients 1 and 2 are the solutions of the system [ where represents the correlation coefficients between the O 3 fields: The ensemble loop starts (see Figure 9): each iteration of this loop produces a new member of the ensemble. 6. The forecast sequence for a single member of the NWP ensemble is retrieved. This sequence will be referred to as "M fields". They have a spatial resolution of 1 km and temporal resolution of 10 min.
7. Growth and decay factors are computed based on the NWP sequence (see section 3.3). 8. The weights for the nowcast-NWP merging are computed (see section 3.4). 9. A sequence of spatially correlated noise images to be used as a shock term in the autoregressive equation are generated. The spatially correlated images A s (t = n ) represent zero-mean, white-noise processes in time (Priestley, 1981;Wei, 2006), correlated in space and having variance Sequences of spatially correlated noise are produced employing a Fourier filter scheme (Pegram and Clothier, 2001a;Bowler et al., 2006;Berenguer et al., 2011;Seed et al., 2013). Extensions of this scheme, such as localized noise characterized by different anisotropy in different locations of the domain, can also be used .
Next, the autoregressive scheme is run to produce a spatiotemporally continuous sequence of images as outlined in Figure 9. The method is iterative and the idea is to use the two most recent available images to produce the next image in the sequence. The new image then takes the position of the most recent image in the sequence and the scheme repeats. The loop is initialized by the two most recent O 3 fields and its first output is the first nowcasting image. The sequence of fields R 0 is produced according to the following AR(2) autoregressive scheme: where n ∈ [1 ∶ N] and represents the displacement of the field estimated by NowTrack. This scheme incorporates an advection step in each iteration, which is expressed by Equation 5a. The two most recent O 3 10-min accumulated radar observations (which are already in the Lagrangian frame of reference of the most recent observation according to step 4) are advected one step ahead using the computed motion field . The advected frames are then placed in Equation 5b to generate the next image in the sequence R 0 ). The new last two images take the position of the most recent images of the sequence (Equation 5c), and the process continues until the desired maximum lead time T max has been reached. At the end of this process, the sequence is exponentiated to make all values positive. This completes the first part of the nowcasting sequence. It is important to understand that both fields P ′ 1 and P ′ 2 are always "centred" in the same frame before Equation 5b is applied to produce the next image. The output of this step is a sequence of N = 36 images. These fields are the basis on which the localization adjustment will be applied later to produce the final nowcast.
At this point, the "localization" loop (see Figure 9) starts. This loop operates on the R 0 fields, applying a localization adjustment, which is now described in detail.

A regular grid of points
is constructed which covers the existing radar domain (observe that the first and last boxes are of smaller size, but it is important to include them for the bilinear interpolations that take place later in the process). Hereafter, this will be called the "coarse grid". The grid length of the coarse grid was chosen to be 32 km in both x and y dimensions. The choice of grid length was decided on the grounds of output performance and computational cost. Our experience suggests that grid lengths in the range 16-64 km function satisfactorily. Consider a single grid point (x di , y di ) of the coarse grid and a box of size 65 × 65 km 2 centred at this point. Notice that boxes at adjacent grid points overlap. We will consider the points on the coarse grid as destination points. The values of two variables are computed for this box: (a) the mean value of the pixels of the box (hereafter image mean flux (IMF) after Pegram and Clothier, 2001a), and (b) the wetted area ratio (hereafter WAR), that is, the fraction of pixels within the box with a value larger than zero. Employing the NowTrack motion field , we can locate the origin coordinates (x oi , y oi ) (frame at time t − ) of destination coordinates (x di , y di ) (frame at time t) and compute the corresponding IMF and WAR at the origin. (If the origin is located out of the domain, IMF and WAR are assigned a zero value). Using the concept of Lagrangian persistence, it is assumed that the values of IMF and WAR are conserved between the origin box (time t − ) and the destination box (time t). At the end of this process, values of IMF and WAR have been assigned at every grid point of the coarse grid and updated at each t. Notice that the first origin frame used is the most recent observation. Mathematically, the operation is expressed as follows: where the formulation IMF[R 0 (x d , t)] represents the IMF of the box centred at location x d of the R 0 image of time t, and the prime in R ′ 0 (x d − , t − ) represents the modification the R 0 field has been subject to under the same transformation in the previous iteration as described later in steps 12 and 13. The outcome of this step is therefore the two grids IMF i and WAR i with i ∈ [1, 32, 64, … ]. 11. The IMFs and WARs are multiplied by the growth and decay factors (see section 3.3). 12. The image R 0 (x, t) produced in step 9 through the autoregressive process (Equations (5a,5b,5c)) now has to get locally adjusted to be consistent with the values of IMF d and WAR d on the coarse grid computed in step 10. A shift-and-scale process, similar to that described in Clothier and Pegram (2002), is used. In particular, for each grid box centred at x di of the coarse grid, the values of two new variables and are computed, so that if (a) is subtracted by all the pixels of the box, (b) all negative values are zeroed, and (c) the outcome is multiplied by , the resulting box of pixels has the desired IMF and WAR (that is, IMF=IMF d and WAR=WAR d ). This way, a value of and is assigned at each grid point of the coarse grid. The outcome of this step is therefore the two grids i and i with i ∈ [1, 32, 64, … ]. 13. Then the shift-and-scale parameters, and , are interpolated bilinearly to a fine grid with grid length equal to 1 km. These maps are defined as K and Λ and have resolution 1 km. Subtracting K from the autoregressed image R 0 , zeroing all negative values of the image, and then multiplying the result by Λ will produce the adjusted image. This process will be referred to as "shift-and-scale" transformation T 1 : R 0 → T 1 R 1 . Observe that the adjusted image has the same spatial resolution as the radar observations, that is, 1 km, while the box size of the coarse grid refers to the scale (64 km) where Lagrangian persistence has been applied. This produces a transformed sequence of images that will be referred to as R 1 . 14. A quantile transformation (probability matching) is applied to R 1 to match the distribution of the resulting image in the sequence with the distribution of the last observation. This, in essence, corresponds to the inverse transform of the normal quantile transform performed in step 4. This produces a transformed sequence of images referred to as R 2 . This transformation affects all the pixels in the image, including those pixels where growth and decay is more intense, but it does not affect the positions where this growth/decay takes place. In essence, despite the fact that the distribution of the most recent observation stays intact, the growth and decay mechanism can reduce or augment precipitation in certain boxes, assigning the precipitation in a way that is consistent with the local orography.
Steps 13 and 14 also help to correct any bias generated, because the zeroes are set to the minimum value of the transformed distribution in step 4. 15. The localization loop repeats to reach the maximum lead time T max . Notice that the iterative adjustment on the 65-km scale, which is described in step 10, causes a propagation of properties (IMF and WAR) in space and time (therefore it operates in a Lagrangian framework). This operation takes place iteratively and results in a gradual advection of the IMF and WAR measures through the terrain, based on the principle of Lagrangian persistence and the underlying motion field. Observe also that the Lagrangian persistence is tweaked to a certain degree in step 11, based on the precomputed NWP-based growth and decay factors (step 7), which operate in a Eulerian framework and are described in detail in section 3.3. 16. The nowcast is merged with the NWP outcome. This produces the final sequence of images referred to as R 3 . The process is described in section 3.4. 17. The process repeats for each member of the ensemble.
A NowPrecip sequence of frames taken every 30 min can be seen in Figure 10. Qualitatively, our experience with NowPrecip can be summarized as follows. (a) The output sequences are very realistic in both spatial variability and temporal continuity, and the lifetimes of patterns are reproduced visually well. They resemble typical observation sequences. (b) There are no visual jumps from the most recent observation and the first nowcasting image at lead time t = +10 min. (c) Merging with the COSMO output takes place seamlessly without obvious discontinuities. (d) The time of production is around 1.5 min and is therefore consistent with the needs of an operational nowcasting system.
A nowcasting ensemble can be produced by repeating this process for a number of members. This can take place either linearly or in a parallel fashion. An example of some members of an ensemble and the corresponding observations can be seen in Figure 11.

Growth and decay
The localized approach of NowPrecip facilitates the estimation of the growth and decay of precipitation. Growth and decay estimation is especially useful for countries with strong orographic gradients, like Switzerland, where interactions between orography and air masses are important.
Research on growth and decay using the most recent observations has been presented before ( The observation is presented in the last panel for comparison. The dashed-line square represents the region that was used for verification research on growth and decay based on multiyear radar archives has been published recently . The NWP precipitation sequences are employed to inform us about the growth and decay in play. For every box of the coarse grid centred at coordinates (x di , y di ), the IMF and WAR at origin (x oi , y oi ) of the NWP frame at time t − and destination (x di , y di ) of the NWP frame at time t are computed. At the end of this process, the values of IMF and WAR at origin and destination coordinates for each grid box and for each of the available 10-min NWP frames in the range between 0 hr and the maximum lead time (in our case +6hr) have been implemented. Therefore the couples (IMF di (t),WAR di (t)) and the corresponding (IMF oi (t − ),WAR oi (t − )) with i ∈ [1, 32, 64, … ] are available.
We have explored two methods to employ this information. Both methods operate in a Eulerian framework.
1. Dynamic growth and decay. For every point of the coarse grid and every lead time t, the following ratios are computed in order to establish local factors of growth and decay: .
[Correction added on 8 April 2020, after first online publication. The typesetting of Equation 7b has been corrected in this version.] These factors are dynamical, that is, they evolve with time at each grid point of the coarse grid. (Observe that, for points where the origin is located out of the domain, both the IMF and WAR growth and decay factors are taken as unity.) Moreover, the sequences of these factors are different between ensemble members, that is, there are different (although typically similar) evolutions of growth and decay. 2. Static growth and decay. First, linear regressions are executed for the variables (IMF di (t),IMF oi (t − )) and (WAR di (t),WAR oi (t − )) for each grid node i separately (a total of N = 36 pairs for each grid point). It is assumed that the intercepts are zero, therefore the two linear regressions are forced to pass from the origin. As a result, the linear regressions provide us with two slopes I i and W i for each grid node of the coarse grid, one for IMF and one for WAR, which represent the average growth and decay factor of that box and do not change with time. There is a second step that may be used to refine this outcome. One may measure the percentages at each box that represent growth and which decay. If growth is far more persistent than decay or vice versa, then this is a box of interest, because it is characterized by persistent growth or decay. One can choose to use the computed rates only for these points, while everywhere else is unity (no growth or decay). If the second step is skipped the method still functions, because at points of no persistent growth and decay the average growth and decay is small. Images of the growth and decay maps can be produced by applying a bilinear interpolation using the IMF and WAR at the grid points.
The net outcome of this process is IMF and WAR growth and decay factors for each grid point of the coarse grid. In the case of dynamic growth, these factors evolve with time, while in the case of static growth and decay they stay unaltered and instead represent an averaged behaviour.
The growth and decay factors enter the autoregressed sequence as follows. During the course of the nowcast simulation, the IMF and WAR for each box and each time frame are produced as usual, but they are now multiplied at each step by the computed growth and decay factor, and before the final bilinear interpolation takes place (step 11 in section 3.2). This way, the statistical properties IMF and WAR of precipitation are adjusted according to the NWP-based growth and decay analysis. In this sense, the concept of seamless forecasting is enhanced, because the nowcasting part can benefit from the NWP not only in terms of a final merging but also in a more decisive way through the evolution. This contribution of the NWP-based growth and decay should have a positive impact on the merging of the nowcast with the NWP product, which happens last in the NowPrecip chain. Now, observe that the orographic effects on the precipitation are captured in part by the motion field, a fact that has been observed before (Li et al., 1995). For instance, motion vectors may get progressively smaller just before locations where orography is associated with decay, and this decrease actually acts as a decay mechanism itself. However, the estimated motion field does not focus on growth and decay per se. The role of a dedicated growth and decay mechanism is essentially to correct the forecast precipitation, especially at locations where growth and decay are systematic, for example, where orography plays a role. The size of the growth and decay factors may be small or large, depending on the extent to which the motion field has already captured the existing growth and decay.
To make this clearer, it helps to have a look at the continuity equation of precipitation (Rinehart, 1981), which is given as The left part of Equation 8 describes the total change of precipitation. The first term on the right describes the local rate of change and this reflects the cloud physical processes that cause generation and dissipation of hydrometeors. The second term on the right, which is called the advective term, represents the advected part of precipitation that has arrived through the motion of clouds by the wind field. It would be elegant for a nowcasting scheme to separate motion from growth and decay efficiently and to treat them separately. This is evidently nontrivial, because of their inherent coupling. Both contributions coexist and interact. The difficulty lies in the fact that in the radar image we only see the superposition of the two effects. The persistence assumption basically expects the second term on the right to dominate significantly over the first. When growth and decay are not very strong, the persistence assumption applies properly and the estimated motion field is adequate. Since, however, the estimated motion field is a combination of both advection and growth and decay, when strong growth and decay take place, the estimated motion vectors end up representing not only the air motion but also, to one extent or another, the underlying growth and decay.
The optical flow scheme described earlier generates the scale that is important for each flow automatically and produces the motion field. Then it analyses local growth and decay to use it as a complement to the motion field. Conceptually, it attempts to describe the two components represented by the two terms in the right part of Equation 8, but since this is never entirely possible, it tries to produce two components that, added together, describe the left part of this equation as accurately as possible. This is particularly beneficial at places where growth and decay happen systematically, for example, at locations characterized by complex orographic terrain.
The computation of growth and decay relies on NWP ensemble members and assumes that the NWP motion is fairly consistent with the motion field produced by our system. This is not a bad assumption in general, although it is not perfectly accurate and it definitely leads to certain errors. More complex schemes than the one described here are feasible. It is possible, for example, to get growth and decay from the most recent observations separately from the growth and decay from the NWP and merge the two seamlessly. Figure 12 shows the maps of growth and decay factors as defined in Equation 7, derived using the static method mentioned above (as described in the technique, the computation relies on information from the available NWP sequence for lead times 0-6 hr). These factors are unitless and show the amount by which the IMF and WAR in each box have to be multiplied to account for the expected growth and decay, as mentioned in step 11 in section 3.2. Factors bigger than 1 denote growth, while factors below 1 denote decay. There is an upper limit for allowable growth and this has been set up to be equal to 3. One can see that growth and decay are, in general, weak at locations not associated with orographic forcing. This is because growth and decay happen everywhere, but are not systematic everywhere. They move with the Lagrangian framework of the precipitation, and consequently average down to small values that reflect only a general trend. However, there are locations associated with the underlying topography, where growth and dissipation are systematic. It is exactly at these locations that the proposed modelling makes a difference. The corresponding sequence of maps for dynamic growth and decay can be seen in Figure 13. In this map, one can see that the growth and decay in the north are not systematic and change with time, while the decay around the coordinates (750,150) is systematic and persists for some hours, disappearing slowly.
The effect of the growth and decay mechanism can be seen in an example. Figure 14 shows a comparison of nowcasts with the growth and decay mechanism on and off. The motion is northwards, while the growth and decay map, which is precomputed (step 7 in section 3.2 and Figure 9) using the static method mentioned above and therefore does not evolve in time, is applied to the sequence of autoregression images, frame after frame (step 11 in section 3.2 and Figure 9). The nowcast without a growth and decay mechanism simply advects the precipitation northwards and makes it pass the Alps and advance towards the Swiss plateau. As one can easily see in the observation, this is not what actually happened. The Alps affect the precipitation evolution and rainfall usually decays, resulting in no precipitation on the Swiss plateau. The growth and decay mechanism analyzes the NWP sequences and finds that there is a strong growth region (in orange) in the southeast and a large decay region (in blue), which essentially delineates the Alps. Applying this growth and decay to our nowcast, in the way described earlier, makes a very strong difference in the outcome. In particular, very small amounts of precipitation can be seen in the northern part of Switzerland. This is because the large static region of decay affects the sequence of autoregression images, decreasing progressively the amount of precipitation they carry as they pass over this area. The growth region is also beneficial in enhancing precipitation in the southeastern part of the map.

Merging with the NWP
There is a broad range of literature on merging techniques of various complexity (Golding, 1998;Atencia et al., 2010;Hwang et al., 2015;Jasper-Tönnies et al., 2018). For instance, Atencia et al. (2010) used a technique where blending is constructed by decomposing the model into new rainfall areas and advection ones, while Hwang et al. (2015) discuss a technique to apply different weights to blend extrapolation forecasts, which are based on intensities and forecast times. A straightforward technique is described here for merging the outcome of nowcasting with that of the NWP. COSMO fields are transformed into a Cartesian grid, cropped to span the exact region of the radar domain and aggregated to 10-min increments. In order to assess the NWP skill, we compute the mean Pearson correlation between the most recent six 10-min radar observations

F I G U R E 12
Maps of growth and decay factors for the variables IMF and WAR using the static scheme described in the text. These maps are bilinear interpolations of the IMF and WAR values used at the coarse grid nodes. These maps are shown here to provide readers with a qualitative idea of the underlying spatial structure of growth and decay. NowPrecip does not use these maps as presented here, instead it uses the values of IMF and WAR at the coarse grid points before it produces each image of the sequence, as described in the text and the corresponding COSMO images. This correlation can be used as a measure of how fast the merging should take place. Intuitively, the smaller , the slower the merging should take place. On the other hand, when NWP is characterized by high skill it is beneficial for the nowcast to merge with it faster. In this sense, instead of performing the merging in a strictly linear way, it is allowed to depend on the skill of NWP outcome. The merged output is a weighted sum of nowcast and COSMO: where R represents the final merged product (sequence R 2 from 16 in section 3.2), F represents the nowcast before merging (sequence R 3 from step 14 in section 3.2), and M represents the COSMO product. The weights are computed using the following scheme: where N = T max ∕ = 36 is the number of frames of the nowcast sequence produced where T max = 360 min is the maximum lead time, and the factor h = 10 −4 is added in case the minimum or maximum is equal to zero. The choice of values for (Figure 10) is empirical. Merging is not applied in the first 10 min in order to keep the skill high; it only begins at a lead time equal to 20 min. Figure 15 shows the process of merging for = 0.0, = 0.2, and = 0.5. For > 0.2, the merging is faster (bigger confidence in COSMO) and for < 0.2 the merging is slower (smaller confidence in COSMO).
The proposed scheme is found to be robust for merging a nowcast smoothly with an NWP outcome. Nowcast ensemble member sequences are merged randomly one-by-one with the corresponding NWP ensemble member sequences. The algorithm is straightforward to code, fast in execution, and can be applied to single-or multiple-member systems without difficulty. Observe that we rely entirely on the correlation of NWP with the observations and not the corresponding one between previous nowcasts and the observations, assuming that the first relates to the second. The main limitation of using the NWP skill in real time is that we do not know the actual shape of the future weighting function, which must be represented a priori. For instance, the ensemble spread is one source of information to help in deriving it. Despite that, numerous simulations using the current scheme suggest that the merging result is qualitatively seamless and without temporal discontinuities, although some smoothing effects can be observed.

F I G U R E 13
Maps of growth decay factors for the variables IMF using the dynamic scheme described in the text. These maps are bilinear interpolations of the IMF values used at the coarse grid nodes. These maps are shown here to provide readers with a qualitative idea of the underlying spatial structure of growth and decay. NowPrecip does not use these maps as presented here, instead it uses the values of IMF (and the corresponding WAR) at the coarse grid points before it produces each image of the sequence, as described in the text decay. The precipitation proceeds from south to north and surpasses the Alps. Top right panel: nowcast with growth and decay. The precipitation, informed by the growth and decay mechanism, does not surpass the Alps. The growth/decay algorithm also motivates growth in the southeastern part of the map. Bottom left panel: the corresponding observation. Bottom right panel: the growth and decay map for IMF. This map is precomputed and static (it was computed for the time stamp 08:30:00, 2018-04-29, and stays unaltered throughout the entire nowcasted evolution (+6 hr) of this run) and applies progressively to the sequence of autoregression images as described in the text (steps 7 and 11 of section 3.2) and in Figure 9 4 VERIFICATION

Verification of NowTrack
Standard techniques (Baker et al., 2011) and two well-known measures are employed for verification of NowTrack: (a) the flow endpoint error (EPE), which represents the Cartesian distance between the true and estimated values (Otte and Nagel, 1994;Zhu et al., 2017), and (b) the angular error (AE), which represents the angle between the true and estimated values (Fleet and Jepson, 1990;Barron et al., 1994). The endpoint error for a pixel at location x, y is defined as where v x and v y represent the true velocity components in the horizontal and vertical directions andv x andv y are the corresponding estimates. Red curve: = 0.0. Black curve: = 0.2. Blue curve: = 0.5. When = 0.2, the merging takes place linearly. When > 0.2, the merging takes place faster than linearly (bigger confidence in COSMO), while when < 0.2 the merging takes place slower than linearly (smaller confidence in COSMO) The angular error is defined as Both synthetic and real fields are used to quantify the accuracy of the geostatistical optical flow scheme. Simple synthetic fields are constructed using the rotation equations: v x = r cos( + x ) + a x , v y = r sin( + y ) + a y .
This scheme simplifies to a circular rotation field for a x = a y = x = y = 0. The parameter denotes the angular velocity. The speeds get larger with distance from the centre of the rotation (see Figure 16, top left panel). In the example used, = 0.04. Experimentation with different choices of the parameter values can produce interesting motion fields. The verification has three steps: (a) produce the synthetic field, (b) apply the synthetic motion field to advect the precipitation field one step ahead, that is, from time t 0 to time t 1 = t 0 + , and (c) use the observed field in t 0 and the synthetically advected field in t 1 as input for the geostatistical motion field estimator to estimate the velocity field. Ideally, the estimated motion field should be the same as the synthetic field used as the truth. The two measures EPE and AE are then used to evaluate the errors.
A similar method that does not employ synthetic fields is also used. The motion field is estimated by NowTrack using frames at times t 0 and t 1 . This motion field is taken as true and is employed to advect the precipitation field from time t 1 to time t 2 . Then one employs the fields at times t 1 and t 2 as input for the optical flow algorithm. The motion field that advects the precipitation field from t 1 to time t 2 is considered as true and the precipitation fields at these two times are now used to reproduce the true motion field.
Results from the method that employed the synthetic field of differential rotation can be seen in Figure 16. In the top left panel, the arrows represent the original true motion field. In the top right panel, the estimated motion field, using the aforementioned procedure, can be seen. The agreement is good in qualitative terms, since the differential rotation has been reproduced to a good extent. The largest differences exist in regions far from nonzero precipitation, something that is expected due to lack of information regarding the actual motion in these areas. However, as mentioned before, NowTrack provides a metric of uncertainty (the kriging variance), which informs users regarding the confidence they should have at every location in the motion field; regions far from nonzero precipitation are characterized by large variance. The bottom left panel shows the endpoint error. This error is very small, in the range of 0-2 pixels at nonzero precipitation locations, and it increases with distance from these locations. A similar picture is seen for the angular error, which is depicted in the bottom right panel. It is noted that small errors do not characterize only regions where the true vectors have small size. Regions with sizeable true motion vectors are also reproduced with excellent accuracy.
Performing the same experiment but using the second method (no synthetic fields) produces results that are also good ( Figure 17). Both EPE and AE are very small in general. Observe that the size of the angular error can be a bit misleading occasionally. When the true motion speed is very small, the motion field estimator may indeed generate a vector of very small speed, but very erroneous direction. In practical terms, this does not have a significant effect on the follow-up extrapolation, because the speed is very small, but the measured AE may appear as large.
Similar behaviour was observed for other examples not shown here. The small errors that characterize NowTrack suggest that it may be possible to use it as a robust scheme to derive information for both large-and small-scale features of precipitation. The example presented was chosen intentionally to be mixed precipitation, that is, characterized by both small-and large-scale features. Despite that, the algorithm managed to produce a very small error even where there are convective cells embedded in the I G U R E 16 Motion verification using synthetic differential motion. Top left panel: the synthetic motion field of differential rotation.
Top right panel: the reproduction of the synthetic field using NowTrack. Bottom left panel: endpoint error map (in km). Bottom right panel: angular error map (in degrees). Observe that the errors are mostly in the range 0-1 km for EPE and 0 • -10 • in regions of nonzero precipitation, but they grow, as expected, progressively larger with distance from these regions background stratiform field. Clearly, further investigation is necessary, but these results suggest that NowTrack succeeds simultaneously in advecting large-and small-scale features for at least some tens of minutes, to a certain extent bridging the gap between areal and object nowcasting systems.

Verification of NowPrecip
Traditional verification metrics are necessary to provide an evaluation of the accuracy of a forecasting system, but additional measures are needed to describe to what extent a nowcast is successful (e.g., Gofa et al., 2018).
An obvious example is when a small-scale feature in the forecast appears displaced by a few kilometres from the position of the observation; then pixel-based verification imposes a double penalty (Ebert, 2008) and the higher the resolution, the worse the situation. A verification of NowPrecip was realized using a variety of measures, which included classic, deterministic ones (Pearson correlation, probability of detection (POD), false-alarm ratio (FAR), and critical success index (CSI)) and measures for neighbourhood verification, like the fractions skill score (FSS) and probabilistic measures (Brier score, ROC curve, reliability diagram, rank histogram). The measures employed for verification are well known, but for completeness they will be defined here briefly. Comprehensive information I G U R E 17 Same as Figure 16, for a different example where now the field that is taken as a reference is realistic, since it has been produced by NowTrack, as explained in the text can be found in standard references (Wilks, 2011;Jolliffe and Stephenson, 2012).

Correlation coefficient between observation and nowcast
(without subtraction of the mean: Zawadzki, 1973). This is a measure that is traditionally used in the deterministic verification of nowcasting systems: where R(x, y, t) represents the 10-min accumulation of precipitation for lead time t at location x, y andR(x, y, t) represents its corresponding forecast value. 2. Probability of detection (POD) or hit rate (HR), false alarm ratio (FAR), False Alarm Rate (FA), and Critical Success Index (CSI). These measures are used for deterministic categorical verification and rely on contingency tables (Stanski et al., 1989), which are constructed by assigning to each grid point of the observations and the forecasts the characterization "event" or "no event", depending on whether or not precipitation has exceeded a certain threshold. The threshold used here is 1 mm⋅hr −1 . These measures are defined as follows: where a is the number of hits (both forecast and observation are classified as events), b the number of misses (observation is an event and forecast a "no event"), c the number of false alarms (observation is a "no event" while forecast is an event), and d the number of correct negatives (both observation and forecast are "no events"). 3. Fractions Skill Score (FSS). The fractions skill score (Roberts, 2008;Roberts and Lean, 2008) is a neighbourhood verification skill score. It has a range between 0 (complete mismatch) and 1 (perfect) and relies on comparisons between boxes of pixels instead of single pixels. As the size of the box used to compute the FSS becomes larger, the FSS improves. The following is the definition. The mean-square error of the observed and forecast fractions for a neighbourhood of length n is where O (n) represent observation fractions and M (n) forecast fractions of points that have exceeded a given threshold within squares of length n, and N x , N y are the number of columns and rows in the domain. The FSS is essentially the MSE score relative to a low-skill reference forecast and is defined as where the reference used for each neighbouring length n is 4. Brier Score (Brier, 1950). This measure is used commonly for probabilistic verification to evaluate ensembles. The Brier score is based on contingency tables; when it is equal to 0, it represents a perfect forecast, and forecasts are worse the larger the Brier score is.
where p(x, y, t + ) is 1 if the event occurred and 0 if it did not, f t (x, y, t + ) = 1 n m ∑ n m j=1p t,j (x, y, t + ), and p t,j (x, y, t + ) is 1 or 0 for events forecast by the jth member of the ensemble at time t for the location (x, y) and lead time ; n is the number of points in the domain, while n m is the number of members of the ensemble.

Relative Operating Characteristics (ROC) curve.
The deterministic ROC score is used to test the discrimination ability of categorical forecasts. Within this context, events are defined as binary, event or no event. Therefore, the ROC score is based on contingency tables, which provide the HR and FA for either deterministic or probabilistic forecasts. The probabilistic ROC scores are calculated for different probability thresholds to form the ROC curve. This curve passes through points (0,0) and (1,1). Forecasts with no skill are indicated by the diagonal line (where HR = FA), while the further the curve is towards the upper left-hand corner, where HR = 1 and FA = 0, the better. The area under the ROC curve (AUC) is actually a measure of separability between classes and the higher it is, the better the system is at predicting zeroes as zeroes and ones as ones. 6. Reliability diagram. The reliability diagram (Hartmann et al., 2002) is a graph of the observed frequency of an event plotted against the forecast probability of an event. It is constructed in similar fashion to the ROC score, but instead of plotting the HR against the FA, the HR is calculated only from the sets of forecasts for each probability separately, and is plotted against the corresponding forecast probabilities. Therefore, in reliability diagrams, comparisons are made against the diagonal, which indicates a perfect score. It is important also to plot the distribution of the forecast probability to show the associated sharpness. 7. Rank histogram. The rank histogram is used to test the ensemble reliability (Anderson, 1996;Hamill and Colucci, 1998;Talagrand et al., 1999;Jolliffe and Stephenson, 2012). It is based on the analysis of how the observed outcomes rank with respect to the corresponding ensemble members. This information provides insight into the average dispersion characteristics of the forecasts. For a reliable ensemble, the members are indistinguishable from observations. Therefore, a reliable ensemble prediction system yields a flat rank histogram. When rank histograms deviate significantly from flatness, the forecasts are characterized by deficiencies.
Now we pass to the verification results. Overall, the verification process produced rather expected results. The outcomes by different measures complement each other and provide a comprehensive picture of the performance of NowPrecip. All deterministic scores have been computed individually for each member of each hour of each event and then these scores were averaged for all available data, that is, all 21 members for each hour of each of the 16 precipitation events (Panziera et al., 2011;Mandapaka et al., 2012;Foresti et al., 2015). Therefore the sample consists of 3,297 members in total. The mean value and bars that indicate the 25 and 75% quantiles of each score's   Figure 18. The corresponding scores for COSMO are also presented as benchmarks.
In general, the scores of each nowcast alone depend on the event itself, the skill of COSMO used for merging, and the time-distance of the hour from the production time of the COSMO output. In general, lower scores characterize nowcasts (a) of convective events, (b) relying on NWP outcomes of low skill, and (c) relying on NWP outcomes that have been produced a long time before the beginning of the nowcast. The verification process took place inside a box smaller than the radar domain. This box extends between the Swiss coordinates (280,-70) km and (910,390) km and includes the entire Swiss area plus large regions of the neighbouring countries. This region can be seen in the bottom right panel of Figure 11. Nowcasts were generated for the whole hours of events presented in Table 1. Each nowcast consists of 21 members. Each nowcast member was merged with a different member of the NWP product. Operationally, the COSMO ensemble is produced twice a day, at 0000 and 1200 UTC, and it becomes available, on average, 85 min later. In the NowPrecip simulations, this lag was taken into consideration. Therefore, nowcasts at time 1300 UTC were merged with COSMO outputs produced 13 hr earlier, at 0000 UTC, since in an operational mode the COSMO output at time 1200 UTC would not have been available yet. The motion fields were produced using 200 points and dynamic growth and decay were used in all simulations.
The results of deterministic skill scores are depicted in Figure 18. In each graph, three curves are depicted: blue represents the NowPrecip output before merging, black represents the NowPrecip output after merging, and red represents the COSMO product. For NowPrecip with merging and COSMO, the bars of the 25 and 75% quantiles of the distributions are also depicted for each lead time. The graph of the correlation coefficient suggests that, on average, there is skill better than that of COSMO for a period up to about 4 hr. POD, FAR, and CSI agree that, on average, there is skill for up to 4 hr. The FSS improves, as expected, with its smoothing length. In particular, the FSS for scale 11 km suggests that there is skill for up to 2 hr. The FSS for 51 km shows that this skill approaches 5 hr. Moreover, there is improvement in the first 2 hr regarding the rate of decrease of FSS 51 km compared with that of FSS 11 km. Intuitively, this may relate to the fact that 51 km is close to the spacing of the grid points of our coarse grid (32 km) and the size of the boxes we have used (65 km). The gain of the merging with COSMO can be seen in all graphs, but mainly in the FSS. Observe that the transition from NowPrecip to COSMO happens seamlessly and there is a quasi-triangular region between the black, red, and blue curves, which gets larger with increasing scale (it is larger for for FSS 51 km than FSS 11 km). This suggests that the gain from the merging process increases with scale. Observe also that, in all the scores, the size of the spread of COSMO skill is relatively constant even at time zero, in contrast with the bars of NowPrecip, which start small, since for early times all members are relatively close to the observation and close to each other. Moreover, the improvement percentages of each of the deterministic scores relative to the COSMO output are presented in Figure 19, employing the formula where S n represents the score of the nowcast, S ref is the score of the reference, which in this case is the COSMO output, and S best is the best value for this score (for instance S best = 1 for POD, while S best = 0 for FAR). More specifically, Figure 19 adds insight to our conclusions and shows the significance of the improvement without the case-to-case variability. The probabilistic scores are presented in Figure 20. The mean Brier score over all the hours of the events saturates in 2 hr, and suggests that there is skill for up to 4 hr compared with COSMO. The reliability diagram, ROC curve, and rank histogram are demonstrated for lead time +1 hr and are compared with the corresponding scores of the COSMO output. The reliability diagram and the ROC curve both suggest a clear improvement over the outcome of COSMO. The curve of NowPrecip in the reliability diagram is close to the diagonal, but still has room for improvement. The sharpness diagram for both NowPrecip and COSMO is also provided for comparison, which demonstrates that NowPrecip probabilities are sharper. In the ROC curve, it is rewarding that FAR is very small for all classes and the NowPrecip curve is pulled towards the upper left corner. This increases the AUC, showing an improved separability between classes over COSMO.
The rank histogram was computed after zeroing the input data below 0.1 mm⋅hr −1 . Moreover, data where both the observation and the maximum value of the 21 forecasts are less than this threshold are not included in the computation. When the observation is zero and the number of zeroes of the forecast members is larger than zero, we randomly assign the rank among the zero forecast positions following the method of Hamill and Colucci (1998). The rank histogram of NowPrecip is very similar to that of COSMO and shows a peak of about 7% at the first bin on the left and a peak of about 12% on the last bin on the very right, therefore underdispersion is observed. These peaks are caused by the inability of nowcasting systems to predict generation from zero or decay to zero at locations where there is no indication in real time that such evolution is about to happen. This effect is intensified for convective cases like the ones mainly employed in the verification process, while the effect is abated in events easier to predict, for example, large-scale stratiform cases. It is expected that the results will improve with an increasing number of ensemble members.

CONCLUSIONS
A new nowcasting system named NowPrecip, currently pre-operational at MeteoSwiss, was presented. It is based on three pillars: (1) the motion field is computed using a geostatistical scheme, (2) there is localized infrastructure, and (3) a growth and decay mechanism is employed.
1. The presented geostatistical motion field generator, NowTrack, has the following advantages.
• The scales of interest of the motion field are computed on the fly. The vectors for a small number of points are computed independently from each other and, after application of quality control to exclude outliers, ordinary kriging is applied to generate the motion field. very intricate localized motions, for instance convection or rotation, as long as they are at the dominant scales obtained by the semivariogram. When a mix of scales is present, the most important one naturally dominates in the estimation of the semivariogram parameters.
• The outlier control excludes unrealistic vectors automatically. The quality control excludes measurements that are inconsistent with the average behaviour of their neighbours. Such discrepancies may be generated, for example, by the presence of very localized growth or decay, which occasionally confuses the motion field estimation process, by producing vectors that represent occasional small-scale convection but do not reflect the bulk motion of precipitation that transports these convective elements on a larger scale.
• Smooth, continuous, and artefact-free motion fields occur in locations where there is no precipitation. This is a characteristic outcome of a kriging-based interpolation and it is critical, because the most recent observed precipitation eventually moves into these regions. Moreover, the motion field is actually probabilistic, since it is accompanied by the kriging variance, although the potential of this feature has not been explored in this work.
2. The platform that underlies the extrapolation is localized and is based on a grid. The points on this grid are 32 km apart. The concept relies on estimating the IMF and WAR at each box, employing Lagrangian persistence in combination with growth and decay information, and producing the nowcasting sequences. This localized platform has the following advantages.

Reliability Diagram
Forecast probability Observed Relative frequency 0 5 10 F I G U R E 20 Probabilistic skill scores. Top left: Brier Score. Top right: ROC curve for lead time +1 hr. Bottom left: Rank histogram for lead time +1 hr. Bottom right: Reliability diagram for lead time +1 hr. The subpanel provides the sharpness distributions. Black represents the NowPrecip outcome after merging, blue NowPrecip before merging, and red represents the COSMO product • Localized behaviour. The concept of localization using a coarse grid is critical for domains characterized by complex terrains, like the one in Switzerland, because in such terrains there are always cases with distinct precipitation patterns in different parts of the domain. This scheme ensures that the statistics of the precipitation in some part of the domain will not affect the nowcast at distant locations.
• The platform is computationally elegant and efficient. The localized treatment introduces a technical basis for treating only a small number of grid points. This makes it easy to incorporate growth and decay into the computation without excessive computational cost.
3. The assessment of growth and decay in nowcasting schemes is a challenge. The proposed scheme is one of the salient features of NowPrecip and utilizes the COSMO output to produce growth and decay information that complements the motion field. The idea is consistent with the concept of seamless forecasting. The NWP outcome relies on far more information than nowcasting, in terms of measurements and theoretical foundation, so it is far more sophisticated concerning the role of numerous variables in the evolution. This information is commonly underutilized by nowcasting systems, but the algorithms described are capable of analysing it and incorporating it into nowcasts. This is particularly useful at locations where the growth and decay is persistent, like regions characterized by complex orography. A verification using 16 multihour events, most of them involving convection, suggests that the method has skill better than COSMO for lead times in the forecast horizon of 0-4 hr. Probabilistically, the verification suggests that the spread of the ensemble members captures reality to a certain extent, although some underdispersion is observed.

FUTURE WORK
NowPrecip has been built with the concept of localization in mind. This arises from needs present in regions with complex terrain like Switzerland, but NowPrecip can equally well be employed in regions with flatter topography. The methods described consist of a stable first version of NowPrecip. There are several points that can be improved or enhanced.
1. The variance of NowTrack can possibly be used to generate perturbations of the motion field. It is conceivable that this variance relates to the errors of the motion field, although this is a claim that has to be investigated and verified properly. To the extent such a relation exists, it can provide ground for including motion perturbations within the scheme. 2. A more complex algorithm may replace the current merging process. A localized version of the presented merging algorithm is currently in testing. According to this, the same quantities as above are computed, but now for each box of our grid individually. This way, merging is performed with different rates according to how well the NWP correlates with the latest radar observations at that location. In locations where NWP performs well, the merging takes place faster than in areas where the performance is poor. Moreover, a complex scheme that attempts to optimize merging, taking into account the spread of the NWP ensemble and nowcasting ensemble, is currently being investigated in MeteoSwiss . 3. Machine-learning capabilities have already been implemented within NowPrecip, but not tested entirely. The localized architecture of NowPrecip makes it naturally machine-learning ready, unlatching the potential for using extrapolation schemes more complex than persistence. In particular, NowPrecip is currently capable of executing extrapolation using on-the-fly a well-known machine-learning technique, decision trees (Hastie et al., 2009;Gareth et al., 2013). It computes IMF and WAR values at each coarse grid point of the most recent observations, collects the data, and accordingly trains a decision tree, which then is used to predict IMF and WAR at the grid points of the coarse grid. This is essentially a nonlinear regression scheme, which has the potential to improve over that of Lagrangian persistence. Observe that it is the small number of grid points of a grid that makes this approach feasible in real time. The same basis can be used for analogue-based machine-learning systems, which rely on multiyear radar archives. Such a system has already been published . 4. Mechanisms that can assist precipitation generation in more targeted ways can potentially improve skill, especially to reduce the number of misses in the rank histogram of Figure 20. This expectation was the main stimulus for developing a method to this end . The scheme in brief is the following. Occasionally, satellite imagery can provide information about locations with relatively high potential for a convective cell initiation, minutes before this gets detected by the radar as a cell. It is useful to generate and incorporate such a forecast cell in a nowcasting output. The method developed generates and evolves convective cells in a realistic fashion. These cells become part of the evolution and interact with their surroundings as integral parts of the evolution. The convective cells are generated stochastically, their evolution can be probabilistic (different for each member of the ensemble), and there is control over their trajectory, size, intensity, anisotropy, and life cycle. An example has been presented in fig. 6 of Sideris et al. (2018). This method proposes a way to use real-time satellite-based prediction on convection initiation in a nowcasting product. The challenge is determining to what extent satellite technology is able to provide solid information on the generation and life cycle of such cells. Currently, some efforts are taking place in MeteoSwiss in this direction. 5. The localized architecture of NowPrecip opens the door to adjusting the precipitation evolution locally based on information coming from other sources. For example, if an independent object-based nowcasting system provides the evolution of a cell's precipitation probability, this piece of information can be inserted quite efficiently in NowPrecip by altering accordingly the IMF and WAR of only those boxes that lie along the expected trajectory of the cell. Therefore, the localized design has the potential to bridge area and object nowcasting systems and present a seamless result within a single platform. 6. The current work aims to present a comprehensive methodology of NowPrecip in a way detailed enough to be reproducible. We chose to accompany this by a multievent overall verification, both deterministic and probabilistic. Additional verification and sensitivity analyses are necessary, especially to investigate all the potential of NowTrack and the growth and decay mechanism, and should be presented in follow-up works. 7. The important concept we attempt to communicate through the growth-decay design is that there is information incorporated within the NWP output that is worth extracting and incorporating inside a nowcast, before the typical final-step merging, and that should be able to assist the seamless paradigm. Incorporation of information derived by the NWP output in a nowcast (before merging) essentially conditions the nowcast for smoother merging. This concept is worth investigating and verifying further. particular to Daniel Leuenberger for providing us with archival data for the completion of this project. We are also grateful to scientists from several institutions, who shared with us their advice, intuition, and experience in numerous workshops and conferences. NowPrecip has been inspired mainly by SBMCast, MAPLE, the String of Beads model, and STEPS, but has used ideas from numerous other works. We acknowledge these contributions and express our gratitude. Loris Foresti and Daniele Nerini were supported by the Swiss National Science Foundation Ambizione project "Precipitation attractor from radar and satellite data archives and implications for seamless very short-term forecasting" (PZ00P2_161316). Special thanks to Isztar Zawadzki and Alan Seed for their expert advice and support. ORCID I. V. Sideris https://orcid.org/0000-0001-7454-7285 U. Germann https://orcid.org/0000-0002-8539-7080