## Abstract

Particle filters (PFs) are Monte Carlo data assimilation techniques that operate with no parametric assumptions for prior and posterior errors. A data assimilation method introduced recently, called the local PF, approximates the PF solution within neighborhoods of observations, thus allowing for its use in high-dimensional systems. The current study explores the potential of the local PF for atmospheric data assimilation through cloud-permitting numerical experiments performed for an idealized squall line. Using only 100 ensemble members, experiments using the local PF to assimilate simulated radar measurements demonstrate that the method provides accurate analyses at a cost comparable to ensemble filters currently used in weather models. Comparisons between the local PF and an ensemble Kalman filter demonstrate benefits of the local PF for producing probabilistic analyses of non-Gaussian variables, such as hydrometeor mixing ratios. The local PF also provides more accurate forecasts than the ensemble Kalman filter, despite yielding higher posterior root-mean-square errors. A major advantage of the local PF comes from its ability to produce more physically consistent posterior members than the ensemble Kalman filter, which leads to fewer spurious model adjustments during forecasts. This manuscript presents the first successful application of the local PF in a weather prediction model and discusses implications for real applications where nonlinear measurement operators and nonlinear model processes limit the effectiveness of current Gaussian data assimilation techniques.

## 1. Introduction

Reliable probabilistic predictions of mesoscale weather events, such as severe convective storms, require accurate estimates of model forecast uncertainty. A large portion of this uncertainty comes from errors in initial conditions, which can grow rapidly over a short period. Forecast errors are typically heterogeneous and depend greatly on the underlying flow, which motivates the use of Monte Carlo data assimilation and forecasting methods for quantifying this uncertainty (Evensen 1994; Houtekamer and Mitchell 2001). For convective-scale weather, which is the focus of the current study, ensemble Kalman filters have proven to be robust methods for assimilating observations and generating probabilistic forecasts (e.g., Snyder and Zhang 2003; Dowell et al. 2004; Caya et al. 2005; Tong and Xue 2005; Aksoy et al. 2009; Dowell et al. 2011). The success of these methods has led to their present use in sophisticated ensemble prediction systems, where they are applied to assimilate routinely collected in situ and remotely sensed radar and satellite observations into high-resolution mesoscale models (e.g., Schwartz et al. 2015; Wheatley et al. 2015; Jones et al. 2016).

Despite recent progress obtained using ensemble-based data assimilation methods, short range forecasts for potentially hazardous mesoscale weather events remain below the level of accuracy needed to provide timely warnings to the public. National agencies, such as the National Oceanic and Atmospheric Administration (NOAA), have dedicated resources toward projects targeting this problem. One example is the NOAA “Warn-on-Forecast” initiative, which involves the collaboration of scientists from a variety of subdisciplines in atmospheric science to achieve operational forecast improvements at time scales ranging from minutes to hours (Stensrud et al. 2009, 2013).

A major challenge for quantifying the probabilistic evolution of convective storms lies in the nonlinearity of the underlying system dynamics. Evidence of the true state comes from relatively sparse in situ measurements, leaving the system poorly constrained by observations. The resulting probability distributions for model analyses and forecasts are non-Gaussian because of the nonlinear evolution of errors between observation times. Densely spaced measurements collected from remote sensing platforms, such as satellites and radars, may be crucial for reducing initial condition uncertainty (van Lier-Walqui et al. 2012). The observed quantities, however, often relate nonlinearly to prognostic state variables in weather models, thus limiting the effectiveness of data assimilation techniques used regularly for weather analysis and prediction (Vukicevic and Posselt 2008; Posselt and Bishop 2012). These methods include ensemble Kalman filters and ensemble–variational hybrids, which typically assume the model dynamics and observation operators are linear. They also assume Gaussian probability distributions for prior and observation error distributions, which are inappropriate for variables such as water vapor and other classes of hydrometeors (Posselt et al. 2014).

In addition to initial condition uncertainty, errors from other sources can create issues for generating accurate and reliable probabilistic forecasts if not accounted for properly. For example, systematic errors in forecast models tend to result in biases for observation-space verifications made during successive data assimilation cycles (Romine et al. 2013). As will be discussed in this paper, these biases can be difficult to distinguish from errors introduced during data assimilation. In particular, corrections made to state variables often do not satisfy the equations of the model used to simulate the dynamical system of interest, which leads to unphysical transient adjustments immediately after initializing weather forecasts; these adjustments are typically referred to as the model’s response to “unbalanced initial conditions.” Dynamical imbalances often come from issues related to sampling error in ensemble-based data assimilation methods (Kepert 2009; Greybush et al. 2011), but they can also come from inappropriate multivariate updates made to model variables that relate nonlinearly to each other. Future progress on reducing initial condition error and model error in weather models may be limited by the assumptions of data assimilation methods used regularly for operational forecasting and research.

Motivated by the challenges listed above, this study uses a new data assimilation technique based on the particle filter (PF) to explore potential benefits of nonparametric filters for mesoscale weather prediction. The method, first introduced in Poterjoy (2016, hereafter P16), assimilates observations serially to approximate the bootstrap PF (Gordon et al. 1993) near observations in state space, while retaining the prior solution outside a radius of influence. In regions between each observation location and the specified localization distance, the local PF produces a quasi-nonlinear update by merging particles together. This update provides a source of imbalance during data assimilation, which becomes less severe as the radius of influence is increased with larger ensemble sizes. Because the filter uses a correlation function to reduce the impact observations have on state variables displaced spatially across the model domain—in a manner similar to localization in ensemble Kalman filters (Houtekamer and Mitchell 2001; Hamill and Whitaker 2001)—the method is called the local PF. Poterjoy and Anderson (2016, hereafter PA16) show that the local PF operates efficiently for high-dimensional geophysical systems and provides benefits over Kalman filtering techniques when linear/Gaussian approximations are invalid. This method is one of several types of localized nonlinear data assimilation techniques that have emerged for problems of moderate to large dimension [e.g., Bengtsson et al. 2003; Lei and Bickel 2011; Frei and Künsch 2013; Van Leeuwen et al. 2015 (chapter 2); Robert and Künsch 2017; Penny and Miyoshi 2016; Chustagulprom et al. 2016]; see Rebeschini and van Handel (2015) and Morzfeld et al. (2017) for an informative review of localization in PFs.

The current study applies the local PF in the Weather Research and Forecasting (WRF) Model (Skamarock et al. 2008) to demonstrate that the techniques outlined in P16 are sufficient for reducing the well-known dimensionality constraint of standard particle filtering methods (Bengtsson et al. 2008; Bickel et al. 2008; Snyder et al. 2008, 2015). In addition to providing one of the first successful tests of a PF-based data assimilation scheme in a weather model, this study also examines deficiencies in ensemble Kalman filtering methods typically used for this purpose. Assimilating observations to analyze and predict mesoscale convective systems presents a case that encompasses the challenges listed above for quantifying errors (Posselt 2016). For this reason, numerical simulations of an idealized squall line in the WRF Model are used to compare the local PF and the ensemble adjustment Kalman filter (EAKF; Anderson 2001). Squall lines are governed by fluid motions over a large spectrum of scales, ranging from small-scale updrafts within convective cells to synoptic-scale features such as frontal boundaries. This study focuses primarily on the role of data assimilation in producing probabilistic representations of convective-scale and system-scale features in mesoscale weather systems, which motivates the use of an idealized model framework. The experiments are carried out using the model setup and observing network described in Sobash and Stensrud (2013, hereafter SS13) and Sobash and Wicker (2015, hereafter SW15), which was initially devised to examine the sensitivity of the EAKF to tunable parameters related to covariance localization and inflation.

The manuscript is organized in the following manner. Section 2 introduces the model and experiment setup used for cycling data assimilation experiments with the local PF and EAKF. Section 3 discusses prior and posterior error statistics and demonstrates advantages of using the local PF for updating cloud variables in the WRF Model. Section 4 presents results from ensemble forecasts initialized with local PF and EAKF posterior members. The last section summarizes the study and discusses future applications of the local PF for model development, atmospheric data assimilation, and forecasting.

## 2. Model and data assimilation setup

### a. Model and observation network

The objectives of this study are to explore potential limitations of Gaussian-based data assimilation techniques for convective-scale weather prediction, and examine the potential of the local PF as a practical alternative to methods currently being used for this application. For these reasons, we perform numerical experiments using a framework with no systematic biases in the model and observing system. Because these experiments use a perfect model, biases introduced from linear/Gaussian assumptions made during data assimilation can be detected more easily than in more complex cases with an imperfect model.

Our experiments adopt the same weather forecast model, truth simulation, and observation network as SS13 and SW15. The model domain is a 396 km 396 km 16 km grid with a 3-km horizontal grid spacing (133 133 horizontal grid points) and 40 vertical levels, most of which are located in the lower troposphere. The area covered by the domain is large enough to contain the entire squall line over a 3-h integration period as it propagates eastward over the course of the simulation. The experiments use version 3.2.1 of the Advanced Research core of the WRF Model with open boundary conditions, WRF single-moment 6-class microphysics scheme (WSM6; Hong and Lim 2006), and a 1.5-order subgrid-scale turbulent kinetic energy scheme. The model contains no terrain, no radiation, and no surface and boundary layer physics.

Initial and boundary conditions for the truth simulation contain 20 m s^{−1} 0–5-km westerly shear and a homogeneous thermodynamic environment taken as the median of 28 proximity soundings from strong mesoscale convective systems as described in Coniglio et al. (2004) (see Fig. 1 of SS13). Figure 1 shows column maximum reflectivity calculated at four times during this simulation. Our truth simulation starts with a line of five discrete cells (Fig. 1a), which are triggered using -K warm bubbles with random perturbations added along the edges to initiate three-dimensional motion. These cells split and merge multiple times before organizing into a mature squall line after 90 min. The environmental conditions allow for deep upright lifting at the leading edge of the gust front, thus allowing the squall line to persist until the end of the 3-h simulation.

As described in SS13, radar radial velocity and reflectivity observations are generated from the truth simulation every 5 min using a radar located at the center of the domain with 14 scan elevations between 0.5° and 19.5°. Errors in these measurements are assigned by adding unbiased Gaussian noise with standard deviations of 2 m s^{−1} and 2 dB*Z* for wind and reflectivity observations, respectively. We then apply the local PF and EAKF with 100 ensemble members to assimilate the observations over the entire life cycle of the squall line. As in SS13, ensemble members are initialized at the beginning of each experiment using the reference profile with random uncorrelated Gaussian noise added to the horizontal wind fields, and keeping the thermodynamic profile the same for all members. These profiles also determine the boundary conditions for each member, which are fixed over the course of the simulation period to maintain a constant source of environmental uncertainty.

### b. Data assimilation methods

The nonparametric filter used in this study is the local PF introduced in P16 and PA16. The local PF operates in a manner similar to the bootstrap PF, except it processes observations serially to generate ensemble members conditioned on observations in the neighborhood of model grid points. This process involves merging sampled and prior members in a manner that considers the posterior mean and variance, given by a matrix of posterior weights, then applying a probability mapping step for higher-order corrections. P16 provides a detailed description of the filter algorithm, along with a schematic illustrating each step. The localized particle update relies on a bootstrap resampling step for choosing which particles to merge with prior particles for forming new posterior particles, which can provide a source of bias for small ensembles. Similar updates are feasible using deterministic approaches as well. For example, Reich (2013) solves an optimal transportation problem to transform prior particles into posterior particles. Nevertheless, the low computational cost of the P16 filter motivates us to use this method in the current high-dimensional application.

The local PF has been included in the NCAR Data Assimilation Research Testbed (DART; Anderson et al. 2009) software suite to take advantage of preexisting infrastructure developed for large geophysical models, such as the WRF Model used in the current study. DART also provides a means of comparing the local PF with ensemble filters applied regularly in atmospheric science, including the EAKF. The DART EAKF is a robust data assimilation technique that is used across multiple disciplines in geoscience (Anderson et al. 2009). In particular, several past studies have used the DART EAKF for convective-scale data assimilation research applications similar to the case adopted for our experiments (e.g., SS13; SW15; Wheatley et al. 2015; Jones et al. 2016). Therefore, this filter provides a good benchmark for testing new data assimilation techniques, such as the local PF applied in this study.

Both the local PF and EAKF require heuristic localization and variance inflation methods to compensate for sampling errors and other systematic errors introduced by approximating the model pdf using a relatively small sample of model states. During the serial assimilation of observations, localization is performed in the EAKF framework by multiplying linear regression coefficients with a Gaspari and Cohn (1999) correlation function, which reduces the state update to zero for grid points greater than a specified distance away from observations. The local PF uses the same empirical correlation function for localization, except it is applied in the calculation of vector weights used during the update [see Eq. (8) of P16]. Both filters also contain parameters that control the amount of ensemble variance, which must be tuned for optimal performance. As discussed in section 2c of PA16, the local PF contains a parameter *α*, which enforces a minimum value for particle weights calculated during data assimilation. This parameter acts to stabilize the filter in situations where localization alone is insufficient for preventing weight collapse (i.e., in cases where observations land outside the span of the ensemble). The *α* parameter effectively reduces the impact of observations that will cause weights to collapse onto a single particle. The EAKF, however, uses the adaptive prior multiplicative inflation technique described in Anderson (2007).

Before arriving at the configurations presented in this study, we performed several experiments to tune the local PF and EAKF parameters. These tests examined the impact of using horizontal localization half widths of 8.91, 10.19, 11.46, and 12.73 km (cutoff 0.0014, 0.0016, 0.0018, and 0.0020 rad in DART), and determined 11.46 km to be the optimal value for both methods. Here, the optimal configuration is taken as the set of parameters that produce the smallest posterior root-mean-square errors (RMSEs). For applications where prior error statistics are close to Gaussian, P16 and PA16 found the local PF and ensemble Kalman filters to exhibit a similar sensitivity to localization, which is consistent with the current study. As in SW15, we keep the vertical localization radius fixed at 3 km; tests examining the sensitivity of vertical localization with the horizontal localization fixed at 11.46 km yielded little sensitivity to this parameter.

The sensitivity of the local PF and EAKF to filter stability/variance parameters was examined in conjunction with the localization sensitivity experiments. The adaptive inflation method used for the EAKF requires users to specify a minimum prior inflation standard deviation and a damping parameter that controls how much the inflation is relaxed to unity each cycle. We chose a prior inflation standard deviation of 0.8 for our experiments, after finding the EAKF to be only mildly sensitive to this parameter. Altering the damping parameter, however, yielded a larger range of results. We tested values between 10% and 50% before arriving at a value of 20%. Likewise, the coefficient *α* in the local PF was tested for values between 0.6 and 0.99. While an *α* value of 0.99 is sufficient for preventing filter divergence, we found to provide the best results.

We also tested both data assimilation methods with and without the additive variance inflation technique introduced in Dowell and Wicker (2009), which adds spatially correlated noise to regions near high-reflectivity measurements. A previous study by SW15 examined the sensitivity of EAKF data assimilation to additive inflation for the same squall-line case used here (see section 2b of their paper for details). Similar to SW15, we use additive inflation with smoothing length scales of 4 km horizontally and 2 km vertically, and error amplitudes of 1 m s^{−1} for wind and 1 K for temperature and dewpoint for data assimilation cycles between 0 and 20 min in the experiment. The added noise helps initiate convection in ensemble members, which start out as laminar flow fields at the beginning of the experiment. After 20 min of cycling data assimilation, we find the local PF to perform best when the error amplitudes are reduced to a quarter of their original values. SW15 used a similar configuration in the first set of cycling EAKF experiments presented in their study. Likewise, the EAKF performs best when additive inflation is replaced completely with Anderson (2007) adaptive multiplicative inflation after this period. The configuration of the EAKF used in the current study deviates from SS13 and SW15, which relied entirely on additive inflation throughout experiments.

## 3. Cycling data assimilation results

This section presents prior and posterior error statistics resulting from 180 min of cycling data assimilation using the local PF and EAKF. To summarize our results, we use volume-averaged quantities calculated within the vicinity of the developing squall line. We first interpolate the model output to equally spaced vertical levels from 0 to 12 km every 500 m, so the verification can be performed using the same uniform vertical coordinate system for each state vector. The vertical levels used for this interpolation are also used as the coordinates for defining the vertical component of the verification volume. To choose the horizontal coordinates, we calculate the column maximum reflectivity of the truth, local PF, and EAKF posterior mean states, then select model grid points where reflectivity values exceed 10 dB*Z* by any three of the model states. We use this volume for verification because it encompasses a portion of the model domain capturing the development and decay of convective cells within the true storm, while also considering regions where spurious convection occurs in local PF and EAKF members. The verification is performed for zonal, meridional, and vertical winds (*u*, *υ*, and *w*), temperature *T*, and mixing ratios of water vapor , cloud water , cloud ice , rain , snow , and graupel .

### a. Prior and posterior mean errors

First, we evaluate the prior and posterior ensemble mean and variance resulting from data assimilation experiments performed over the life cycle of the squall line. Figure 2 shows time series plots of RMSEs and spread, calculated at each 5-min data assimilation cycle. Prior and posterior values are plotted together at each cycle time, thus producing the sawtooth pattern in the figure. We omit the first 15 min of the experiments as both methods rely heavily on additive inflation to initiate convective cells in each ensemble member.

In general, both methods produce reasonably accurate analyses over the course of the experiment; for example, most cycles contain posterior RMSEs for *u* and *υ* that are less than the observation errors for radar winds. The resulting mean errors are also lower than the 50-member EAKF experiments presented in SS13 when verified over the same volume (not shown), which is expected for the 100-member experiments performed here. The comparison reveals that the EAKF produces lower RMSEs than the local PF for all variables except (Fig. 2). One reason for this result is the quasi-linear evolution of cells between the data assimilation cycles presented here, which helps maintain linear relationships between the observed and unobserved quantities in the prior, thus making the assumptions of the EAKF appropriate.

Despite producing lower RMSEs than the local PF, the EAKF tends to overestimate error variances in several state variables; namely, *u* and *υ*. This overestimation is partly due to the environmental profiles used to specify model boundary conditions, which are sampled from an unbiased error distribution and held constant throughout the experiments. The mean of this sample contains only a small amount of bias, which comes from sampling error in the 100-member ensemble. This relatively unbiased sample leads to the transport of ensemble variance into the storm region with no corresponding increase in mean RMSEs, thus contributing to a slight overestimation of errors by both the local PF and EAKF ensembles for *u* and *υ*. The larger overestimation of variance by the EAKF is due mostly to suboptimal uncertainty estimates that come from assimilating reflectivity observations. Though not shown, omitting reflectivity observations from the EAKF experiment yields ratios of RMSE to ensemble spread that more closely resemble that of the local PF experiment.

A previous study by P16 noted that the relative performance of the local PF compared to ensemble Kalman filters depends on both the nonlinearity of the application and the ensemble size used to perform the data assimilation. For the current application, we hypothesize that the local PF would eventually produce smaller RMSEs than the EAKF with increasingly larger ensemble sizes, which follows from the low-dimensional model experiments in P16. Exploring this limit for convective-scale data assimilation in the WRF Model, however, is not the goal of the current study. We also suspect the 100-member ensemble used here is likely too small to capture high-order moments of the posterior error distribution accurately, as noted in Miyoshi et al. (2014). Nevertheless, a closer evaluation of the local PF and EAKF experiments reveals benefits of the nonparametric filter for the application presented here. We will show in section 4 that the lower posterior errors of the EAKF ensemble do not yield lower forecast errors than the local PF, because of physical imbalances that occur when using a linear filter for a nonlinear/non-Gaussian application.

### b. Probabilistic cloud analysis

The application in this study presents a case where variables associated with cloud features may exhibit large deviations from Gaussianity. For non-Gaussian applications, verifying the mean and variance alone is insufficient for quantifying filter performance. To illustrate a typical situation encountered when applying ensemble filters to generate cloud analyses, we present examples using the local PF and EAKF to assimilate a single radar reflectivity measurement. For this demonstration, we choose an observation located near a high-reflectivity region within a convective cell 90 min into the simulation. Figure 3 shows horizontal and vertical model cross sections through the observation to be assimilated (location indicated by the red “×” in the figure); reflectivity and storm-relative winds are plotted on each cross section to provide a three-dimensional overview of the convective cell. The height of the observation is close to the freezing level in the truth simulation (about 3.8 km in altitude), where prior ensembles contain a mix of frozen and unfrozen hydrometeors. For the single-moment microphysics scheme used in this study, the prior representation of clouds in the model comes from ensemble predictions of rain, snow, graupel, cloud ice, and cloud water mixing ratios.

First, we examine the prior ensemble representation of joint pdfs for reflectivity and hydrometeor mixing ratios (blue markers in Fig. 4). The prior ensemble comes from 5-min forecasts run from posterior local PF members 90 min into the cycling data assimilation experiments. Because mixing ratios are bounded at zero, the marginal probabilities of each quantity used in this example are guaranteed to be non-Gaussian. The relationships between hydrometeors and reflectivity are also nonlinear, which further limits the effectiveness of Gaussian approximations for these joint pdfs.

We create an observation (indicated by the “y” on the *x* axis in Fig. 4) by adding a 2-dB*Z* error to the truth (yellow markers in Fig. 4), then assimilate this observation using the local PF and EAKF. The resulting posterior samples (red markers in Fig. 4) demonstrate limitations in the Gaussian approximation used by the EAKF (left column of Fig. 4). The most obvious deficiency of the EAKF is the large number of posterior realizations found in regions far from where the prior indicates nonzero probability. For example, most posterior EAKF members lie outside the narrow region of high joint probability for reflectivity and (Fig. 4c), because of the linear correction applied to prior members. The EAKF also generates many negative values for posterior , , and , which are physical impossibilities.^{1} In practice, negative mixing ratios generated during data assimilation are set to zero before model integration, but the correction is made without considering the dependence of other variables, such as temperature and moisture. The members resulting from the local PF assimilation of this measurement resemble the bootstrap PF solution, which generates a posterior ensemble by sampling from the prior with replacement based on weights assigned to each member (Gordon et al. 1993). Provided that forecast members are samples from the true prior error distribution, the bootstrap PF will converge to the Bayesian posterior solution as ensemble size goes to infinity.

Next, we look at updates made by the two filters to variables separated geographically from the observation. For this experiment, we again examine mixing ratios located at the model grid point indicated in Fig. 3, but use an observation displaced 2.5 km below these values. As in the previous example, samples from the joint pdf of 5-min forecast reflectivity and hydrometeor mixing ratios (blue markers in Figs. 6 and 7) indicate non-Gaussian distributions for prior errors. Figures 5 and 6 also show posterior EAKF and local PF members with and without localization, respectively, to demonstrate its effect on posterior updates. Localization acts to create a posterior ensemble comprising members that share properties of both the prior and unlocalized posterior members, which can reduce the physical fidelity of the resulting updates (Kepert 2009; Greybush et al. 2011). For the EAKF, localization acts to reduce the linear correction coming from assimilating the reflectivity measurement (comparing the first columns in Figs. 5 and 6). The impact of localization in the local PF framework is less trivial, because weights assigned to members during data assimilation are localized instead of regression coefficients (see P16). Without localization, weights far from the observation are identical to the weights at the location of the observation. With localization, weights are tapered to their prior value at grid points far from observations, so it is possible for members with effectively no weight before localization to have a much larger weight afterward. This process can cause some members to shift back to regions of high prior probability after localization (cf. the second columns in Figs. 5 and 6). For example, localization causes posterior values of to move from 0 g kg^{−1} to a secondary region of high prior probability between and g kg^{−1}. The localized updates maintain some of the nonlinear characteristics of the bootstrap PF within the localization region, which helps provide a more realistic depiction of cloud features within the squall line for analysis members.

Nonlinear physical relationships between model state variables present another major challenge for linear/Gaussian data assimilation approaches. The process of generating posterior samples that satisfy the equations of the dynamical model often requires sampling from a non-Gaussian multivariate probability density. Figure 7 shows examples using prior and posterior samples from the joint pdf of *T* and hydrometeor mixing ratios; both sets of variables in each example are again located at the “×” in Fig. 3. The samples shown in this figure demonstrate the complex dynamics represented by the WRF microphysics scheme. For example, WSM6 melts cloud ice as soon as *T* is greater than 0°C, thus producing zero values for about half the members. Likewise, is large for members with *T* > 0°C, but remains close to zero for all prior members above 0°C because of heterogeneous freezing. Mixing ratios for snow and graupel show a slightly more gradual transition between temperature regimes, especially in the *T* > 0°C regime, where frozen particles melt over a deep layer in the model. Because of the non-Gaussian joint pdfs, posterior EAKF members generated by assimilating the reflectivity measurement in Fig. 3 often result in model states that are highly improbable or impossible. In addition to producing many members with negative mixing ratio values, several posterior EAKF members contain large mixing ratios of frozen quantities in the *T* > 0°C regime, and large in the *T* < 0°C regime. Integrating the WRF Model from these states would inevitably require adjustments within cloud regions early in the forecast, which can lead to rapid changes in the concentration of hydrometeors. The resulting adjustments can produce sudden increases or decreases in diabatic heating rates, which feed back into vertical accelerations induced by convective cells in the storm. The posterior local PF members, however, retain prior relationships between mixing ratios and *T* at the location of the measurement (right column of Fig. 7), which helps alleviate some of the issues raised in this demonstration.

Reflectivity fields produced by posterior members provide one means of examining the cumulative effects of using the EAKF and local PF to update cloud variables over many successive data assimilation cycles. Figure 8 shows maximum column reflectivity for the first 10 members of the EAKF posterior (top two rows) and local PF posterior (bottom two rows) from the last data assimilation cycle at 180 min; values are plotted in increments of 5 dB*Z* between 5 and 65 dB*Z*, and the true 5-dB*Z* contour is plotted for reference. The largest qualitative differences between local PF and EAKF members exist in regions where significant deviations from Gaussianity are likely to occur—in this case, along the trailing edge of the squall line where prior ensembles contain a mix of relatively high mixing ratios and zero mixing ratios. EAKF members tend to contain spurious reflectivity on the periphery of convective cells during development, which follows from the suboptimal update of hydrometeor mixing ratios within this region. Likewise, the local PF members contain far fewer regions of spurious convection near the developing storm and yield more realistic depictions of the squall line than the EAKF members.

The positive impact of the local PF over the EAKF, summarized subjectively in Fig. 8, is quantified using rank histograms in Fig. 9. Rank histograms indicate the frequency with which the true model state occurs in bins formed by ranking members in ascending order (Anderson 1996; Hamill and Colucci 1996; Talagrand et al. 1997). Any deviations from a uniform histogram indicate possible deficiencies in the probabilistic representation of variables by an ensemble. For this verification, we populate the histograms using variables located within the leading edge of the squall line at 15-min intervals. The temporal frequency of variables helps reduce correlations between samples, which are assumed to be independent. Though not shown, the histogram shapes also do not change significantly when grid points are subsampled from different parts of the squall line, and at different times. This step ensures biases from various parts of the squall line do not cancel to provide a false impression of good probabilistic skill. The verification is performed for the same variables and altitude used for the demonstration introduced in Fig. 3, which contains nonzero values for each hydrometeor variable at various times and locations during the simulation. As expected, the histograms reveal noticeable errors in the probabilistic representation of mixing ratios by the EAKF ensemble. For example, and contain substantial biases, which come from the suboptimal update of these variables using a linear filter. On the other hand, the local PF produces relatively flat histograms for the variables shown here, thus demonstrating its ability to produce reliable probabilistic representations of non-Gaussian quantities during cycling data assimilation.

## 4. Forecast results

In this section, we examine the skill of forecasts initialized from posterior local PF and EAKF members. We perform this comparison on ensemble forecasts run every 15 min between the 60- and 120-min cycling times. These initialization times are representative of major changes in the life cycle of the squall line, after mean posterior RMSEs become near stationary. The 15-min spacing between initialization times helps reduce correlations between the set of forecast errors used for this verification, which provides six sets of forecasts with approximately independent errors. We run each set of ensemble forecasts for 60 min, which is long enough to determine which data assimilation technique produces more accurate forecasts for the fast-evolving squall line, but short enough to run several sets of forecasts to the end of the 180-min truth simulation. The 60-min forecast period also allows more than enough time for transient physical adjustments to dissipate after initializing the model with local PF and EAKF members.

We verify the forecasts by calculating biases and RMSEs of the state variables, then average the statistics over the same area used in section 3 for the posterior verification. We then average these values over the six sets of predictions to summarize the forecast performance of the two data assimilation techniques. Figure 10 shows mean 60-min forecast bias and RMSEs for the local PF and EAKF experiments. These values are plotted as a function of altitude from the surface to 12 km to compare vertical profiles of errors through the squall line. Despite producing systematically larger posterior RMSEs than the EAKF experiment (Fig. 3), the local PF ensemble yields smaller 60-min forecast RMSEs and bias for nearly all variables and model levels. The remainder of this section will focus on possible reasons why forecasts initialized with local PF members produce more accurate forecasts than EAKF members over the course of the experiment.

First, we examine how errors grow in time by comparing volume-averaged RMSEs every 5 min during the forecasts (Fig. 11). As Fig. 2 demonstrates, the EAKF initially contains lower posterior errors than the local PF for all variables except . Nevertheless, errors in the EAKF forecasts grow much faster than the local PF errors over the 60-min period. Differences in error growth rates between the two sets of initial conditions tend to vary across model variables. For example, the smallest differences in error growth rate occur for *T* and horizontal wind variables. Though not shown, joint prior error distributions between these variables and the observed radar observations are approximately Gaussian. Therefore, both the local PF and EAKF are capable of producing posterior estimates of *u*, *υ*, and *T* that are relatively close to physically balanced model states. These variables are also governed mostly by larger-scale dynamics in the model, which have smaller error growth rates than variables that are predominantly determined by convective-scale processes. For the remaining variables—namely, , , , and *w*—EAKF forecasts yield much larger error growth rates than the local PF over the first 0–30 min. The rapidly growing errors for microphysics variables and vertical motion in the EAKF case indicate a higher transient response to model initial conditions, which comes from unphysical updates made by the filter. In this case, the larger short-range forecast errors following model initialization are a direct consequence of the Gaussian assumptions made by the EAKF (see section 3b).

In Fig. 12, we compare vertical RMSE profiles for , , diabatic heating rate (in K s^{−1}), and *w* between 0 and 10 min into the forecasts. With the exception of , these quantities undergo rapid error increases early in the forecast, due to the low predictability of cloud processes and the transient adjustments that occur after model initialization. The latter of these error growth mechanisms is stronger for the EAKF case, as the states initialized with the local PF tend to remain closer to the true model state during the forecasts. We include in this comparison because it presents a major distinguishing factor between the two data assimilation methods. Of the quantities verified in this study, is the only variable for which the EAKF consistently produces larger posterior RMSEs than the local PF (Fig. 12a). The elevated errors in the EAKF experiment result mostly from a systematic dry bias at low altitudes and a moist bias near the melting level (about 3.8 km), as shown in Fig. 13. The vertical structure of bias resembles the dewpoint biases shown in Fig. 6 of SW15 for the same squall-line case. SW15 also noted strong *T* biases in their EAKF experiments, but were able to decrease these biases substantially by decreasing the amount of additive inflation. The EAKF in this study does not use additive inflation after the first 20 min of cycling; as a result, no substantial *T* bias exists in the EAKF prior and posterior fields. The EAKF also produces larger biases than the local PF for cloud particles, rain, and frozen hydrometeors; we do not show these biases because they contribute to only a small fraction of the total mean RMSE. Because the local PF does not produce the same large biases in , we suspect this result comes from the cumulative effects of updating hydrometeor mixing ratios improperly with a linear filter. One such process, as described in section 3b, follows from the production of excess frozen hydrometeors below the freezing level during data assimilation, which can translate into higher through sublimation, or melting followed by rapid evaporation. The excess cloud particles quickly change phase between data assimilation cycles, leading to a positive bias with no corresponding large bias in frozen or unfrozen hydrometeor mixing ratios. We examine the spatial distribution of these biases by averaging posterior mean errors in a gust-front-relative framework over the data assimilation cycles used to produce Fig. 13. Biases are calculated at the approximate altitude of the freezing level (Figs. 14a,c) and a level located 2 km below (Figs. 14b,d). This figure shows that biases are more pronounced along the periphery of the squall line, where entrainment of dry air accelerates the evaporation of cloud particles. A large region of negative biases occurs directly below the level of excess , likely because cloud particles are evaporating at higher altitudes before reaching this region. Large biases also appear in regions where data assimilation has little immediate effect on the posterior solution because few or zero observations exist nearby. This result leads us to suspect the biases are likely caused by residual effects of suboptimal cloud variable updates in the leading cells. The local PF yields a bias pattern similar to the EAKF for the two altitudes compared here, except values are much smaller and not as concentrated near the western boundary of the squall line as the EAKF case.

The large positive moisture bias helps accelerate error growth for other hydrometeor variables in the EAKF forecasts. For example, the EAKF produces larger RMSEs than the local PF by 5 min into the integration (Fig. 12b). These errors grow fastest near the 2–4-km layer, where the positive bias is largest. Errors in hydrometeor variables translate into errors in diabatic heating rates, which lead to an amplification of *w* errors (Figs. 12c,d). This process presents one plausible explanation for why the EAKF forecasts have larger *w* RMSEs than the local PF forecasts after 10 min of model integration. The EAKF forecasts for *u*, *υ*, and *T* tend to take about 20 min to yield larger RMSEs than the local PF initialized forecasts. This period is approximately equal to the life cycle of clouds in the forecasts; therefore, we suspect these errors reflect the time required for errors at the convective scale to have a significant impact on the system-scale evolution of the squall line.

Based on our analysis of short-range forecast errors, the 60-min forecast benefits demonstrated by the local PF likely come from more accurate moisture analyses and an improved preservation of physical balances when the prior pdf is non-Gaussian. This hypothesis prompted a second set of data assimilation experiments, where the interval between observations is increased from 5 to 20 min. In the latter case, the model state is less constrained by observations than the configuration used throughout this manuscript, which allows errors to grow more nonlinearly with time. Results from the 20-min cycling experiments are summarized using volume-averaged RMSEs (dashed lines in Fig. 11) and profiles of horizontal mean 60-min forecast RMSEs and bias (Fig. 15). As expected, assimilating observations less frequently exacerbates issues related to non-Gaussianity for linear filters, which further decreases the ability of the EAKF to produce accurate forecasts. The resulting verifications show clearer advantages of using the local PF over the EAKF for predicting the evolution of the squall line.

The 20-min cycling scenario is more comparable to real convective weather applications than the 5-min configuration used earlier in this study. For reference, the National Severe Storms Laboratory (NSSL) real-time prediction system assimilates radar, satellite, and surface observations in 15-min cycles (Wheatley et al. 2015; Jones et al. 2016). The cycling frequency used by the experimental NSSL system is a large improvement over convective-scale prediction systems run operationally over the United States, such as the hourly cycled NOAA High-Resolution Rapid Refresh (Benjamin et al. 2016). Therefore, the limitations of using linear filters for real-time convective-scale weather applications are likely more severe than the example shown here.

## 5. Discussion and conclusions

This study presents initial tests of the local PF for data assimilation in a weather prediction model. The application consists of a developing squall line simulated by the WRF Model, which is observed every 5 min by a dense network of observations. Starting from a random ensemble of WRF Model states, radar velocity and reflectivity measurements are assimilated sequentially over the course of the storm’s life cycle using the local PF and a commonly used EAKF data assimilation technique. Both methods use the NCAR DART software package for performing filter update steps of the data assimilation cycles. Results are verified using prior and posterior model states and 60-min ensemble forecasts run at various times during the experiment. To the best of our knowledge, this study presents one of the first successful applications of a PF-based data assimilation method in a weather model.

Using only 100 members, the local PF provides accurate probabilistic analyses of the squall-line simulation examined in this study. Though not presented in the manuscript, the local PF also provides reasonably accurate posterior estimates using as few as 50 members; experiments using 25 members resulted in filter divergence after several cycles. These experiments demonstrate that the local PF can operate effectively for high-dimensional geophysical applications at a cost comparable to ensemble data assimilation techniques currently used. Because joint observation-model space prior error distributions are close to Gaussian for many of the variables, the 100-member EAKF experiment produces smaller posterior RMSEs than the 100-member local PF experiment for most variables verified in this study. Nevertheless, ensemble forecasts initialized with local PF members are more accurate than EAKF initialized forecasts by 10–20 min. The advantages of the local PF come from an improved physical representation of posterior members and smaller biases in hydrometeor variables during cycling data assimilation. In particular, continuous cycling of the EAKF leads to a large water vapor bias in the trailing stratiform region of the squall line, which comes from unphysical updates of cloud variables. One common example is the large number of negative hydrometeor mixing ratios that occur in posterior members, which must be set to zero before model integration. A more subtle case, discussed in section 3b, is the excess amount of frozen hydrometeors the EAKF produces well below the freezing level; likewise, excess rainwater tends to be found well above the freezing level. Posterior local PF members, however, retain the same properties as prior model states near each measurement assimilated during the sequential processing of observations. Merging particles in the localization region can act to destroy dynamical balance during data assimilation, but we find balances to be better maintained by the local PF than by the EAKF. Therefore, advantages of applying a nonlinear filter to update cloud variables in the squall line are large enough to overcome advantages of the EAKF for the remaining Gaussian-distributed state variables in 5-min cycling data assimilation experiments. The benefits of the local PF become clearer when observations are assimilated every 20 min instead of 5 min, which allows errors to grow more nonlinearly with time.

Results from this study encourage future testing of the local PF as a data assimilation method for real geophysical applications. In particular, the local PF may be beneficial for assimilating measurements that pose challenges for Kalman filters or variational data assimilation strategies because of nonlinearity. Similar to the reflectivity observations used in this study, cloudy radiances present one example where both nonlinearity in the measurement operator and nonlinearity in the model dynamics limit the ability of Gaussian filters typically used to assimilate these observations. In addition, the local PF benefits from its flexibility in specifying a pdf for observation likelihood calculations (Poterjoy and Anderson 2016). The current study uses observations with Gaussian errors, but assigning skewed or bounded distributions for measurements with known non-Gaussian errors is trivial in a PF framework. These factors may have important impacts for the assimilation of remotely sensed data provided by dense satellite and radar observation networks. Developing strategies to assimilate these observations more effectively may be crucial for producing significant gains in numerical weather prediction.

Last, the local PF provides a non-Gaussian data assimilation framework for investigating errors in model parameterization schemes. Gaussian assumptions made by the EAKF in this study often lead to biased short-range forecasts for hydrometeor variables. Forecast errors of this form are difficult to distinguish from errors caused by deficiencies in the microphysics schemes themselves, which complicates strategies for diagnosing errors in cloud models. By reducing biases introduced during data assimilation, the local PF presents a means of isolating sources of systematic bias that may occur when assimilating real observations. The presence of unknown model errors also introduces challenges for applying the local PF for numerical weather prediction, because forecast members may not resemble samples from the true prior error distribution. This problem can be avoided by adopting techniques used by PFs that rely on carefully chosen proposal densities, such as the equivalent-weights PF and the implicit PF (van Leeuwen 2010; Chorin et al. 2010). Incorporating model error in a non-Gaussian framework presents one future area of research for this filter.

Nonparametric filters, such as the local PF, may provide powerful tools for investigating unknown error sources in modeling and observing systems and making better use of remotely sensed data. Future work will explore the potential of the local PF for real data assimilation applications and exploit the advantages of this method for improving parameterization schemes in models.

## Acknowledgments

This research is sponsored by the NCAR Advanced Study Program. Any opinions, findings, and conclusions or recommendations expressed in this publication are those of the authors and do not necessarily reflect the views of the National Science Foundation. This study benefited greatly from discussions with Louis Wicker. We would also like to acknowledge high-performance computing support from Yellowstone (ark:/85065/d7wd3xhc) provided by NCAR’s Computational and Information Systems Laboratory, sponsored by the National Science Foundation.

## REFERENCES

*Probability and Statistics: Essays in Honor of David A. Freedman*, D. Nolan and T. Speed, Eds., Vol. 2, Institute of Mathematical Statistics, 316–334.

*Pushing the Limits of Contemporary Statistics: Contributions in Honor of Jayanta K. Ghosh*, B. Clarke and S. Ghosal, Eds., Vol. 3, Institute of Mathematical Statistics, 318–329.

*13th Conf. on Probability and Statistics in the Atmospheric Sciences*, San Francisco, CA, Amer. Meteor. Soc., 51–56.

*Proc. ECMWF Workshop on Predictability*, Reading, United Kingdom, ECMWF, 1–25.

*Nonlinear Data Assimilation*. Frontiers in Applied Dynamical Systems: Reviews and Tutorials Series, Vol. 2, Springer, 118 pp, doi:.

## Footnotes

^{a}

The National Center for Atmospheric Research is sponsored by the National Science Foundation.

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

^{1}

Approximations made during the update and mapping steps of the local PF can also lead to negative mixing ratios, but this result is far less severe and occurs less frequently than in the EAKF experiment.