## Abstract

Nowcasts (short-term forecasts) of heavy rainfall causing flash floods are highly valuable in densely populated urban areas. In the Collaborative Adaptive Sensing of the Atmosphere (CASA) project, a high-resolution X-band radar network was deployed in the Dallas–Fort Worth (DFW) metroplex. The Dynamic and Adaptive Radar Tracking of Storms (DARTS) method was developed as a part of the CASA nowcasting system. In this method, the advection field is determined in the spectral domain using the discrete Fourier transform. DARTS was recently extended to include a filtering scheme for suppressing small-scale precipitation features that have low predictability. Building on the earlier work, Stochastic DARTS (S-DARTS), a probabilistic extension of DARTS, is developed and tested using the CASA DFW radar network. In this method, the nowcasts are stochastically perturbed in order to simulate uncertainties. Two novel features are introduced in S-DARTS. First, the scale filtering and perturbation based on an autoregressive model are done in the spectral domain in order to achieve high computational efficiency. Second, this methodology is extended to modeling the temporal evolution of the advection field. The performance and forecast skill of S-DARTS are evaluated with different precipitation intensity thresholds and ensemble sizes. It is shown that S-DARTS can produce reliable probabilistic nowcasts in the CASA DFW domain with 250-m spatial resolution up to 45 min for lower precipitation intensities (below 2 mm h^{−1}). For higher intensities (above 5 mm h^{−1}), adequate skill can be obtained up to 15 min.

## 1. Introduction

Flash floods can cause significant property damage and loss of lives. Therefore, nowcasts (short-range forecasts) of severe rainfall causing such events are highly valuable for the society. Providing early warnings is especially important in densely populated urban areas. Weather radars are ideal for this purpose because of their high spatial and temporal resolution, for which the typical ranges in operational radar networks are 1–5 km and 5–15 min. Thus, the use of weather radars has become increasingly important in urban hydrology (Liguori et al. 2012; Thorndahl et al. 2017).

In the Collaborative Adaptive Sensing of the Atmosphere (CASA) project, an experimental X-band radar network has been deployed in the Dallas–Fort Worth (DFW) metroplex with over 7 million people. This area frequently experiences hazardous phenomena such as heavy rainfall, flash floods, hail, and tornadoes. As of 2018, the network consists of eight radars, and the typical spatial and temporal resolutions of the reflectivity and QPE products are 250 m and 1 min, respectively (Chandrasekar et al. 2018).

The Dynamic and Adaptive Radar Tracking of Storms (DARTS; Ruzanski et al. 2011) is an extrapolation-based method that was developed as a part of the CASA nowcasting system. In this method, estimation of the advection field is formulated as a linear least squares problem in the spectral domain using the fast Fourier transform (FFT), and a kernel-based method is used for extrapolating the radar echoes. More recently, DARTS was extended to use a filtering scheme for small-scale precipitation features having low predictability, which led to development of Scale Filtering DARTS (SF-DARTS; Pulkkinen et al. 2019). In this method, the scale filtering is done by decomposing the discrete Fourier transform (DFT) representation of the reflectivity field into radial frequency bands and applying autoregressive AR(2) models in the spectral domain.

The limitation of DARTS and SF-DARTS is that they produce only deterministic nowcasts. However, it was shown in Pulkkinen et al. (2018) that the nowcasts in the CASA DFW domain are subject to large uncertainties even for short lead times (less than 15 min). Estimates of the uncertainties are highly desirable especially when using the nowcasts as inputs to hydrological models (Thorndahl et al. 2017). The sources of uncertainty include growth and decay of precipitation and uncertainty in the advection field estimation as well as measurement errors in the input data (Villarini and Krajewski 2010; Liguori and Rico-Ramirez 2014). The standard approach to modeling these is to add stochastic perturbations to a deterministic nowcast in order to produce a nowcast ensemble, from which probability distributions of future precipitation intensities can be obtained (Hohti et al. 2000; Bowler et al. 2006; Berenguer et al. 2011; Dai et al. 2015).

The feasibility of the Short-Term Ensemble Prediction System (STEPS) methodology (Bowler et al. 2006) originally developed for continental-scale nowcasting was evaluated in the CASA DFW domain in Pulkkinen et al. (2018). While the results were promising and showed the advantages of probabilistic nowcasts, the computational complexity of STEPS was observed to be too high for operational use in the target domain, where nowcasts need to be generated every 1 min at a high spatial resolution. One of the main reasons for this was the use of a computationally intensive variational optical flow method for determining the advection field.

To meet the above requirements, the objective of this study is to develop Stochastic DARTS (S-DARTS), a probabilistic precipitation nowcasting method built on the DARTS/SF-DARTS methodology. The key principle is to formulate the nowcasting model in the spectral domain, which leads to high computational efficiency. This is a fundamental difference compared to a majority of the existing stochastic methods (Bowler et al. 2006; Berenguer et al. 2011) with the only exception known to the authors being the phase stochastic method (PhaSt; Metta et al. 2009). In these methods, FFT is applied to obtaining the power spectrum or scale-based decomposition of the reflectivity field, after which the inverse FFT is applied and the nowcast model is applied in the spatial domain.

In S-DARTS, stochastic perturbations are added to the reflectivity and advection fields in order to model the uncertainties in the deterministic SF-DARTS nowcasts. Thus, S-DARTS can be viewed as a hybrid of SF-DARTS and STEPS. Realistic perturbations are generated by respecting the spatiotemporal structure of the observations. A novel feature of S-DARTS is applying a vector autoregressive VAR(2) model for the temporal evolution of the advection field. Additionally, a postprocessing technique is employed to ensure that the nowcasts have the same statistical properties with the observed reflectivity field. A summary of differences between the DARTS, SF-DARTS, S-DARTS, STEPS, and PhaSt methods is given in Table 1.

This paper is organized as follows. The proposed methodology is described in section 2. Performance evaluation and demonstration of using S-DARTS in the CASA DFW domain is given in section 3. Concluding remarks and discussion are given in section 4. Finally, technical implementation details of S-DARTS are given in the appendix.

## 2. Methodology

A schematic view of the S-DARTS methodology is shown in Fig. 1. The input is a time series of radar-measured reflectivity fields, and the output is a stochastically generated ensemble of nowcasts. Descriptions of the individual components are given in the following subsections.

### a. Advection field estimation and extrapolation

In DARTS (Ruzanski et al. 2011), the advection field estimation is done by numerically solving the equation

where $z\u2061(x,y,t)$ denotes the reflectivity field in space and time and $u\u2061(x,y)$ and $\upsilon \u2061(x,y)$ denote the *x* and *y* components of the advection velocities. Differently to the conventionally used correlation or variational techniques (Laroche and Zawadzki 1994; Li et al. 1995; Bowler et al. 2004), the solution is in DARTS done in the spectral domain. This is done by applying FFT to a time series of observed reflectivity fields *z* of dimension $L\xd7L\xd7T$ and then solving for the DFT representations of the advection field components *u* and *υ*.

For computational reasons, the solution of the advection equation is done up to the desired DFT truncation number denoted by *M*. Additionally, the DFT of the reflectivity field is truncated to the range $Nxy$ and $Nt$ in space and time. Denoting the truncated DFTs of the advection field components by $UM$ and $VM$, a discrete approximation of Eq. (1) in the spectral domain is written as (Ruzanski et al. 2011)

where $kx,ky\u2208\u2061[\u2212Nxy,Nxy]$, $kt\u2208\u2061[\u2212Nt,Nt]$, the operator $\u22c6$ denotes discrete convolution, and

Determining $UM$ and $VM$ then amounts to obtaining a least squares solution of the resulting linear equations. For the details of our implementation of DARTS, see Pulkkinen et al. (2019).

For extrapolating the radar echoes, we use the semi-Lagrangian backward interpolate once (IO) scheme described in Germann and Zawadzki (2002). The method is implemented by using bilinear interpolation and three iterations for each time step. The choice of this method instead of the spectral extrapolation used in SF-DARTS is dictated by the fact that here the DFT coefficients of the forecast reflectivity field are not truncated, which is necessary for an efficient implementation. Instead, the portion of the Fourier spectrum suppressed by the SF-DARTS model is filled with stochastic perturbations (see section 2c).

### b. Scale decomposition

The scale decomposition used in S-DARTS is adapted from Pulkkinen et al. (2019). This is done by applying FFT to the $L\xd7L$ reflectivity field and dividing the Fourier wavenumber range $[0,L/2]$ into *n* radial bands $B1,B2,\u2026,Bn$.

The weights for the Fourier wavenumbers $k=\u2061(kx,ky)$ and bands $Bj$ are given by

where $u|k|,j$ are normalization constants to ensure that for a given radial wavenumber |**k**| the weights sum to one and $g|k|,j$ are Gaussian functions of the form

### c. The forecast model for the reflectivity field

In SF-DARTS (Pulkkinen et al. 2019), autoregressive AR(2) models are applied to frequency bands of observed reflectivity fields in the spectral domain. Since predictability is proportional to spatial scale, small scales having low predictability are progressively suppressed with increasing lead time. The S-DARTS model extends SF-DARTS by adding stochastic perturbations. The purpose of this is to simulate nowcast uncertainties and also to fill the portion of power spectrum that is suppressed by the AR(2) process.

Adapting the approach of Pulkkinen et al. (2019), we apply an AR(2) model to each frequency band $Bj$ and sum the contributions from each band to obtain the final forecast field. Applying the AR(2) models in the spectral domain allows separation of frequency bands for computational efficiency, which cannot be done when the models are applied in the spatial domain (e.g., Bowler et al. 2006).

Separate models are used for the deterministic and stochastic components. The forecast equation is written as a weighted sum of these so that

where

and

for $j=1,2,\u2026,n$ and $(kx,ky)\u2208Bj$ and $\Delta t$ denotes the time interval between the observed reflectivity fields.

Here we consider only half of the Fourier spectrum and restrict $kx$ to the range $[0,\u230aL/2\u230b+1]$, since the DFT of a real-valued input is symmetric (Pardo-Igúzquiza and Chica-Olmo 1993; Sundararajan 2001). A detailed description of estimating the parameters $\psi j,1$ and $\psi j,2$, initialization of the models, and generating the perturbation term $\epsilon Z$ is given in the appendix.

The multiplier $\Psi |k|,\epsilon \u2061(t)$ for the stochastic term defined by Eq. (6b) is chosen as

where

and the operator

denotes the variance of a two-dimensional field *z* contained in the band $Bj$ expressed in terms of power spectral density (von Storch and Zwiers 2003). The above choice of $\Psi |k|,\epsilon $ ensures that the portion of the power spectrum suppressed by the deterministic AR(2) model, Eq. (7), is filled with stochastic perturbations obtained from Eq. (8) so that the resulting forecast field obtained from Eq. (6) has approximately the same power spectrum with the observed one at the nowcast start time $t0$.

The AR(2) model defined by Eq. (7) is applied in Lagrangian coordinates. This is imposed upon the initialization of the model by applying the inverse FFT and extrapolating $z\u2061(\u22c5,\u22c5,t0\u22122\Delta t)$ and $z\u2061(\u22c5,\u22c5,t0\u2212\Delta t)$ to the initial time $t0$. The semi-Lagrangian method described in section 2a is used for this purpose, after which the reflectivity fields are transformed back to the spectral domain.

For each time $t>t0$, the final forecast field $z\u2061(\u22c5,\u22c5,t)$ is obtained from Eq. (6) by applying the inverse FFT and the method described in section 2e to ensure that the output field has the correct statistics. Finally, the semi-Lagrangian method is applied to extrapolate the resulting field from *t*_{0} to *t*. For computational efficiency, the extrapolation is done for each time step incrementally by using the previously computed displacement vectors (i.e., the total motion along the trajectory) from time $t0$ to $t\u2212\Delta t$.

### d. The forecast model for the advection field

A novel feature of S-DARTS is that a multivariate generalization of the above AR(2) model is applied to the advection field. This approach is motivated by our hypothesis that filtering unpredictable small-scale features in the advection field with increasing lead time improves forecast skill. The model also provides a framework for adding stochastic perturbations.

In the proposed approach, a vector autoregressive VAR(2) model is applied to the DFT representation

of the advection field, where $kx,ky\u2208\u2061[\u2212M,M]$ (see section 2a). We assume that the mean advection velocity over the domain remains invariant throughout the forecast, and we subtract the mean value before applying the model. Since *M* is typically small (i.e., less than 10), the computations can be carried out in the spectral domain in a highly efficient manner. For this reason, we also simplify the above VAR(2) model so that the frequency decomposition described in section 2b is not applied.

As in section 2c, separate models are used for the deterministic and stochastic components. The formulas are given by

and

where $\Phi 1$ and $\Phi 2$ are 2 × 2 matrices, $W^\u2061(kx,ky,t)=\u2061(0,0)$ for $k=\u2061(0,0)$, and $\epsilon W$ is a perturbation field of dimension $(2M+1)\xd7\u2061(2M+1)$. The methods for estimating the parameters $\Phi 1$ and $\Phi 2$ and generating the perturbation fields $\epsilon W$ are described in the appendix.

The perturbed advection field is obtained by summing the deterministic and stochastic terms given by Eqs. (13) and (14) and the mean-value term

This yields

which is supplied as input to the semi-Lagrangian extrapolation (see section 2a) after zero padding and the inverse FFT. The purpose of the matrix $\Phi \epsilon \u2061(kx,ky,t)$ is to rescale the perturbations to correspond to the observed error distribution of the VAR(2) model. A detailed description of determining the elements of this matrix is given in the appendix.

### e. Postprocessing

The underlying assumption in S-DARTS is that there is no systematic growth or decay in the precipitation field throughout the forecast. Thus, we require that the forecast reflectivity fields have the same probability distribution with the most recently observed one, and also that perturbations are not added into areas of no precipitation. To satisfy these requirements, we apply a three-step masking, shifting, and quantile mapping procedure to the output of the AR(2) model described in section 2c.

The masking and shifting are done using an extended version of the method described in Seed (2003) for deterministic nowcasts. The superlevel sets are defined as

for the perturbed and unperturbed forecast fields *z* and $z\u02dc$ obtained from Eqs. (6) and (7) via the inverse FFT, respectively. Then, for a given contour value *c*, connected regions $Rc,i$ of the set $M\u02dcc$ are identified by using the algorithm described in Wu et al. (2005). These are included as candidates for precipitation areas if $Rc,i\u2229Mc\u22600$.

The optimal contour value $c*$ is obtained by applying a bisection iteration such that the set

has approximately the same number of elements as the number of pixels in the observed reflectivity field $z0$ exceeding the 20-dB*Z* threshold. The masked and shifted forecast field is then obtained from

which has approximately the same fraction of pixels exceeding the 20-dB*Z* threshold as $z0$.

In effect, the above approach avoids generating perturbations far from areas that are identified as no precipitation according to the deterministic nowcast. Since the cutoff between precipitation and no precipitation areas is based on intensities of the perturbed nowcast, this allows generating realistic nowcast fields with no visible artifacts at boundaries of precipitating areas (see Fig. 10 in section 3). This would not be the case if the precipitation mask were taken directly from the deterministic nowcast as in Seed (2003) or from the observed reflectivity field as in Pulkkinen et al. (2018).

As the final step, the conditional CDF of the forecast field $z\u02dc\u2032$ above 20 dB*Z* is mapped to the conditional CDF of $z0$. That is,

where *F*_{0} and *F* denote the conditional CDFs of $z0$ and $z\u02dc\u2032$, and the superscript −1 denotes inverse function.

## 3. Performance evaluation

In this section we evaluate the forecast skill of S-DARTS in the CASA DFW domain with different ensemble sizes. We also present examples involving the use of probabilistic nowcast products.

### a. Description of the data and the parameters

As of 2018, the CASA radar network deployed in the DFW region shown in Fig. 3 consists of eight dual-polarization X-band radars (Chandrasekar et al. 2018). The reflectivity mosaics used in this study were generated by compositing clutter-filtered and attenuation-corrected reflectivity scans. The composition was done by interpolating the volumes from four elevation angles into a latitude–longitude constant-altitude plan position indicator (CAPPI) grid with spatial resolution of 250 m and temporal resolution of 2 min. The size of the resulting grid is 1121 pixels × 1121 pixels. Additionally, 20 dB*Z* was used as the threshold value between precipitation and no precipitation, and reflectivities below 20 dB*Z* were set to 15 dB*Z* for the nowcast computation.

The reflectivity mosaics were collected during 12 convective events in the CASA DFW domain (Pulkkinen et al. 2018). The total length of the dataset is 42 h (approximately 1260 reflectivity mosaics with the 2-min time resolution). A summary of the selected events is given in Table 2.

The parameters used in the experiments are listed in Table 3. These values represent a trade-off between nowcast skill and computational performance. Computation times of DARTS were observed to be highly dependent on the number of DFT coefficients *N*_{xy}, *N*_{t}, and *M*, but experiments with different values (30–100, 3–6, and 2–15, respectively) did not show significant differences in nowcast skill. Therefore, the relatively small values were chosen for best computational performance. The experiments were done on a Linux system with two six-core Intel Xeon E5-2630 processors running at 2.6 GHz.

Since the observed and forecast precipitation intensities are given as reflectivities (dB*Z*), the corresponding rainfall intensity values are listed in Table 4. These numbers are based on a generic *Z*–*R* relation (Marshall and Palmer 1948), and thus are not optimized for the CASA DFW radar network. Nevertheless, they serve as a rough approximation of the rainfall intensities to ease the interpretation of the reflectivity values.

### b. Verification metrics

The relative operating characteristic (ROC) curve (Jolliffe and Stephenson 2003) measures the ability of a probabilistic forecast to discriminate between precipitation and no precipitation with respect to the given reflectivity threshold. Given a set of $np$ probability thresholds, the ROC curve is constructed by plotting the probability of detection (POD) against the false alarm rate (POFD). For an ideal forecast, the curve passes through the upper-left corner (i.e., 100% hit rate and 0% false alarm rate). The area under the ROC curve can be used as a summary statistic. In the following, we use $np=100$ with evenly distributed probability thresholds between 0 and 1.

The rank histogram (Hamill 2001) measures correspondence between the spread of the ensemble and the observed uncertainty. For each pixel in the forecast grid, the values of the members $1,2,\u2026,n$ are ranked from the lowest to highest, and a pooled histogram with values $h1,h2,\u2026,hn+1$ is constructed by assigning each verifying observation a bin that it falls into among the ensemble members [ties are handled by using the method described in Hamill and Colucci (1998)]. The first and last bin are assigned for observations being below or above all members, respectively. For a forecast ensemble whose distribution is consistent with the observations, the histogram is flat with no observations falling into the first or last bin. From the rank histograms, we compute the percentage of outliers

with different ensemble sizes *n*.

The reliability diagram (Bröcker and Smith 2007) measures correspondence between predicted exceedance probabilities of precipitation intensities and their observed frequencies. For a perfectly reliable nowcast, the curve lies on the diagonal. Typically a reliability diagram is accompanied with a histogram of the sample sizes for each probability interval used in computation of the diagram. The forecast is sharp if a majority of the predicted probabilities are assigned to 0 or 1, meaning that the forecast gives definite yes or no predictions with no uncertainty. Here we divide the probability range [0, 1] into *k* intervals, compute a histogram $hi$, $i=1,2,\u2026,k$, of the predicted probabilities, and define sharpness as the fraction

of values in the rightmost bin and the total number of forecast probabilities (the first bin containing zero probabilities is not included).

The above verification statistics were aggregated over the events listed in Table 2. This was done by computing a nowcast every 10 min such that the start and end times of the nowcasts (46 min with 2-min time steps) were within the time range of the event. For each lead time, the nowcast-verifying observation pairs (*n* nowcast values and one verifying observation for each pixel) were then pooled into the same dataset, from which the statistics were computed. Pixels, where the verifying observations and all ensemble members had zero precipitation intensity, were excluded from the computations.

### c. Verification with different ensemble sizes

The ensemble size plays a major role in the computational requirements of S-DARTS and its ability to produce probabilistic nowcasts. Thus, we conducted a set of experiments to determine the number that gives a reasonable nowcast skill while keeping the computational requirements low. As the verification metrics, we use the area under the ROC curve and the percentage of outliers (OP).

Figure 4 shows ROC areas plotted with respect to lead time with reflectivity thresholds 20 and 40 dB*Z* for different ensemble sizes *n*. Increasing the ensemble size gives an improvement in the ROC area, and the improvement becomes more noticeable with increasing reflectivity threshold and lead time. For 20 dB*Z* and a lead time of 45 min, the relative improvement from 6 to 96 members is 4%, and for 40 dB*Z*, we obtain a 10% improvement. We can also observe that increasing the ensemble size from 48 to 96 members gives a nonnegligible improvement in the ROC area only for the 40-dB*Z* threshold (up to 5% compared to less than 1% with the 20-dB*Z* threshold). This indicates that 48 members is sufficient unless nowcasts of precipitation extremes are desired.

Figure 5 shows OP for S-DARTS nowcasts with different ensemble sizes *n* and lead times. The results show very large variability between different ensemble sizes. The average OP for ensemble size *n* = 6 is approximately 45%, whereas it is reduced to under 15% when *n* is increased to 96, and a further reduction can be expected when *n* is increased over 100. In addition, Fig. 6 shows that excluding the lower and upper ends, the rank histograms corresponding to *n* = 48 are highly uniform for all lead times (10, 20, and 30 min shown here), indicating the ensemble spread is sufficiently large to represent the uncertainty in the nowcasts.

Another observation from Fig. 5 is that for all ensemble sizes, the OP statistic first increases until 4–6 min, but then slowly decreases with respect to lead time. At the beginning of the nowcast, the ensemble spread is too small because of uncertainties not accounted for by the S-DARTS model. This is possibly explained by the fact that we use a rather coarse estimate of the advection field (see the parameter *M* in Table 3), leading to loss of small-scale variability. The relative importance of this source of uncertainty decreases with lead time, as the uncertainties originating from the temporal evolution of precipitation intensities becomes dominant.

### d. Verification with high precipitation intensities

The ability to produce nowcasts of higher precipitation intensities (e.g., over 10 mm h^{−1}), that can cause flooding and other hazardous phenomena, is essential for the practical applicability of S-DARTS in the CASA DFW domain. Therefore we carried out a set of experiments to assess the forecast skill of S-DARTS with different intensity thresholds.

The experiments were done by using ensemble size *n* = 48 that gives a good trade-off between computational performance and nowcast skill. With this setting, the computation time of a 1-h nowcast on our test system with 2-min time steps is approximately 3 min. With *n* = 96 the computation time increases to 6 min, which is no longer computationally feasible in the CASA DFW domain where rapidly updating nowcasts (i.e., less than 5 min) are desired.

Figure 7 shows ROC areas plotted with respect to lead time with different reflectivity thresholds and *n* = 48. The nowcast skill decreases with lead time, and the rate of decrease is highly dependent on the chosen threshold. While the ROC area remains above 0.82 for thresholds up to 25 dB*Z* for lead times up to 45 min, a rapid drop can be observed when the threshold is increased to 40 dB*Z* and above. For such values, after 15 min the ROC falls below 0.72, which can be considered as the limit for nowcasts to be usable (Buizza et al. 1999).

Figure 8 shows reliability diagrams plotted for different lead times with 20- and 35-dB*Z* reflectivity thresholds. For the 20-dB*Z* threshold, the forecast exceedance probabilities show high reliability with reasonably small deviations from the observed frequencies. However, for the 35-dB*Z* threshold, the reliability drops rapidly with respect to lead time, being at a reasonable level only for lead times up to 10 min. The exceedance probabilities are overestimated for all lead times. This is possibly explained by the fact that perturbations are generated uniformly over the whole domain without accounting for local variability in the statistics of the precipitation field. Consequently, high precipitation intensities can be generated into areas where they are unlikely to occur (see Figs. 10 and 11).

As shown in Fig. 9, the loss of reliability with increasing reflectivity threshold is accompanied by loss of sharpness. For thresholds up to 25 dB*Z* with lead times up to 30 min, the sharpness is above 0.25 (i.e., the forecast probabilities of exceeding the threshold are close to 100% over 25% of the time). However, when the threshold is increased to 35 dB*Z*, the sharpness falls below 0.2 only after 10 min.

The results shown in section 3c suggest that a significant improvement to the above results cannot be expected by increasing the ensemble size above 48 members. Nevertheless, Figs. 7–9 show that usable probabilistic nowcasts of higher precipitation intensities (e.g., over 35 dB*Z* corresponding to 5 mm h^{−1}) can be generated up to 15 min, which is sufficient for decision-making and preparatory actions. The above observations do not necessarily indicate a deficiency of the S-DARTS model. In fact, they reflect the low predictability of small cells containing high precipitation intensities, which makes nowcasting such phenomena a very challenging task for any method.

It is worthwhile to note that the results presented for S-DARTS show several improvements compared to the STEPS method used in Pulkkinen et al. (2018) particularly for higher reflectivity thresholds. With 25 STEPS ensemble members and reflectivity threshold of 35 dB*Z*, the ROC area was observed to fall below 0.6 after 30 min in Pulkkinen et al. (2018). As seen from Fig. 4, the ROC area for the corresponding S-DARTS ensemble with 40-dB*Z* threshold remains close to 0.7 after 30 min. In addition, the S-DARTS nowcasts show improved sharpness. For instance, with the 30-dB*Z* threshold, the sharpness of the STEPS nowcasts was observed to fall below 0.05 after 20 min in Pulkkinen et al. (2018), but Fig. 9 shows that the corresponding number for S-DARTS is over 0.2. These improvements are explained by the two main differences between S-DARTS and STEPS: statistical postprocessing of nowcasts and a more advanced perturbation scheme for the advection velocities.

### e. Illustrative examples

Figure 10 shows the initial reflectivity field, the unperturbed SF-DARTS nowcast and two members of an S-DARTS nowcast ensemble during event 10. In this case, the precipitation pattern is moving east along the advection field. We can observe from the ensemble members that adding the Gaussian perturbations to the reflectivity field produces realistic patterns that respect the spatial structure of the observations. In addition, perturbing the advection velocities produces precipitation patterns that are displaced from the locations given by the deterministic nowcast, and the displacement is spatially variable.

The main added value of S-DARTS to the deterministic DARTS method currently used in the CASA nowcasting system is the ability to provide exceedance probabilities for the given precipitation intensities. For a given grid cell, the exceedance probability is defined as the fraction of ensemble members above the given threshold. This is illustrated in Fig. 11 showing exceedance probabilities during event 10 for reflectivity thresholds 20, 30, 35, and 40 dB*Z*.

## 4. Conclusions

A novel combination of two methodologies: SF-DARTS and STEPS was proposed as a probabilistic precipitation nowcasting method for the CASA nowcasting system. Adapting the spectral advection field estimation and autoregressive scale filtering from SF-DARTS and stochastic perturbations from STEPS, the key elements of these two methods were combined. The resulting method, called S-DARTS, implements the autoregressive model and stochastic perturbations in the spectral domain. Another novel feature of S-DARTS is application of a vector autoregressive model for the temporal evolution of the advection field.

Verification experiments were carried out in order to evaluate the performance of S-DARTS in the CASA DFW domain at 250-m spatial resolution. Based on the area under the ROC curve, the reliability diagram and sharpness plots, it was shown that for lower precipitation intensities (below 2 mm h^{−1}), S-DARTS can produce skillful probabilistic nowcasts of convective rainfall up to 45 min. For higher intensities (above 5 mm h^{−1}), adequate skill can be obtained up to 15 min. Another observation is the improved nowcast skill compared to the STEPS method used in Pulkkinen et al. (2018). This is attributed to two main enhancements of S-DARTS compared to STEPS: statistical postprocessing of nowcasts and a more advanced perturbation scheme for the advection velocities.

It was also shown that ensemble sizes between 24 and 48 members give a good ability to represent nowcast uncertainties. However, larger ensembles of up to 100 members may be considered when nowcasts of precipitation extremes are desired. Using 96 ensemble members, the percentage of verifying observations falling outside the ensemble was observed to be below 15% for lead times above 10 min. Together with the highly uniform rank histograms, the conclusion is that the large sample limit of the S-DARTS ensemble is unbiased and is well able to capture the variability in the spatiotemporal evolution of precipitation fields.

S-DARTS is feasible for implementation in the CASA nowcasting system using a standard desktop computer. The initialization phase takes less than 30 s, and computing a single ensemble member of a 1-h nowcast with 2-min time steps takes approximately 1 min on one processor core. Provided that enough processor cores are available, computation of the ensemble members can be parallelized without significant additional computational cost.

The S-DARTS model assumes that the spatial distribution of precipitation intensities remains constant throughout the forecast. However, because of rapid systematic growth and decay patterns in small domains such as the CASA DFW, significant biases can be present in the nowcasts. Therefore, modeling the temporal evolution of parameters of the intensity distribution, as proposed in Berenguer et al. (2011), could improve forecast skill. Another limitation of the model is the assumption of spatial homogeneity of the precipitation field, which can be questioned even for small domains such as the CASA DFW (see Fig. 10). Thus, the forecast skill could be improved by generating the perturbations based on localized statistics, as proposed in Nerini et al. (2017).

## Acknowledgments

This research was supported by the Finnish Academy of Science and Letters via the Foundations’ Post Doc Pool and the U.S. National Science Foundation Hazards SEES program (Grant AGS-1331572).

### APPENDIX

#### Implementation Details of the S-DARTS Model

##### a. The frequency decomposition

The weight functions introduced in section 2b are defined by

$g0,j=0$ for all $j=1,2,\u2026,n$ and

where

the number of levels *n* = 8 and *L*_{0} = 3 (Pulkkinen et al. 2019).

The normalization constants $u|k|,j$ in Eq. (4) are chosen as

##### b. The AR(2) models

The AR(2) model parameters $\psi j,1$ and $\psi j,2$ in Eqs. (7) and (8) are estimated by using three most recently observed reflectivity fields. This is done by using time-lagged autocorrelation coefficients

at the nowcast start time $t0$ for $l=0,1,2$ and the Yule–Walker equations (Pulkkinen et al. 2019; Priestley 1981).

When applying the AR(2) models, the frequency bands $Bj$ are truncated such that

and the weights $w|k|,j$ outside the bands $Bj$ are set to zero (i.e., excluded from the computations).

The initial terms of the deterministic AR(2) models defined by Eq. (7) are taken from the DFT representations of the most recently observed reflectivity fields. That is,

where $j=1,2,\u2026,n$ and $l=0,1$. Initialization of the stochastic models defined by Eq. (8) is done by generating two random fields $\epsilon Z\u2061(kx,ky,t0)$ and $\epsilon Z\u2061(kx,ky,t0\u2212\Delta t)$ using the method described in section d of the appendix and setting

for $j=1,2,\u2026,n$.

##### c. The VAR(2) model

The parameters $\Phi 1$ and $\Phi 2$ of the VAR(2) models defined by Eqs. (13) and (14) are estimated from the generalized Yule–Walker equations (Wei 2006):

where

and $w\u2061(t)=\u2061[u\u2061(t),\upsilon \u2061(t)]$ denotes the spatial domain representations of the advection fields that are computed for $t0\u2212l\Delta t$, with $l=0,1,2$. Similarly to Eq. (A7), the initial terms of the VAR(2) model are obtained from two most recently computed advection fields.

The elements of the scaling matrix $\Phi \epsilon $ in Eq. (16b) are given by

where $\lambda 1$ and $\lambda 2$ denote the eigenvalues of the covariance matrix $\Gamma 0$ in decreasing order, and the columns of the matrix $V$ contain the eigenvectors of $\Gamma 0$ corresponding to these eigenvalues. These are also the eigenvectors of the covariance matrix of $w\epsilon \u2061(t)$, since the perturbations are generated so that they have the same covariance matrix as the initial advection field (see section d of the appendix).

The above scaling functions $\alpha 1\u2061(t)$ and $\alpha 2\u2061(t)$ are chosen as third-order polynomials that depend on the lead time $t\u2212t0$ according to

The coefficients $ci,j$ are obtained by fitting $\alpha i\u2061(t)$ to differences between the deterministic VAR(2) model predictions [obtained from Eq. (16a) after zero padding and the inverse FFT] and the observed advection field along the eigenvectors contained in the matrix $V$. The fitting is done based on data from a time window of $Tw$ hours preceding the start time $t0$.

##### d. Generation of the perturbation fields

The perturbation fields $\epsilon Z$ in Eq. (8c) are generated so that they have the same spatial correlation structure with the observed reflectivity field. Based on the Wiener–Khinchin theorem (Kay and Marple 1981; von Storch and Zwiers 2003), this is done by using the observed power spectrum as a filter for white noise. This approach allows generation of perturbations with anisotropic structures such as rainfall bands (Niemi et al. 2014) that are relevant for larger scales (e.g., over 50 km) in the CASA DFW domain (Ruzanski and Chandrasekar 2011).

The reflectivity perturbations are generated by sampling a $\u230aL/2\u230b\xd7L$ field of Fourier phase angles $\theta i$ independently from the uniform distribution $U\u2061(0,2\pi )$ such that the symmetry condition

is satisfied. This is followed by computation of

for $|k|>0$, where the multiplication by $|Z\u2061(kx,ky,t0)|$ ensures that the perturbations have the same power spectrum with the observed reflectivity field $z\u2061(t0)$. Additionally, we set $\epsilon Z\u2061(0,0)=0$ in order to generate perturbations with zero mean. With condition (A13), the inverse DFT of the resulting perturbation field is real valued and normally distributed (Pardo-Igúzquiza and Chica-Olmo 1993). Thus, our approach is essentially the Fourier integral method (FIM) of Pardo-Igúzquiza and Chica-Olmo (1993), but instead of using a parametric covariance function, this information is obtained directly from the observed power spectrum.

Assuming that the errors of the VAR(2) model applied to the advection field are normally distributed, the perturbation fields $\epsilon W$ in Eq. (14c) are generated by using a generalization of the above method. The components of the perturbation field

are obtained from

for $i=1,2,3,4$, where $kx\u2208\u2061[0,M],ky\u2208\u2061[\u2212M,M]$, $\theta U~U\u2061(0,2\pi )$ such that condition (A13) is satisfied and

With the above definitions, $\epsilon U$ and $\epsilon V$ have the same amplitudes and phase differences with the initial advection field components $U\u2061(\u22c5,\u22c5,t0)$ and $V\u2061(\u22c5,\u22c5,t0)$, which implies the same correlation and cross-covariance structure (von Storch and Zwiers 2003).

## REFERENCES

*Forecast Verification: A Practitioner’s Guide in Atmospheric Science.*John Wiley and Sons, 254 pp.

*Spectral Analysis and Time Series.*Academic Press, 890 pp.

*The Discrete Fourier Transform: Theory, Algorithms and Applications.*World Scientific Publishing, 374 pp.

*Statistical Analysis in Climate Research.*Cambridge University Press, 484 pp.

*Time Series Analysis: Univariate and Multivariate Methods.*2nd ed. Addison Wesley, 614 pp.

## Footnotes

© 2019 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).