JOBNAME: MWR 00#00 2022 PAGE: 1 SESS: 8 OUTPUT: Sat Dec 24 08:35:26 2022 Total No. of Pages: 21/ams/mwr/0/mwrD220108

MONTH 2022

1

KUROSAWA AND POTERJOY

A Statistical Hypothesis Testing Strategy for Adaptively Blending Particle Filters andEnsemble Kalman Filters for Data Assimilation
KENTA KUROSAWA
a

a

a,b

AND JONATHAN POTERJOY

Department of Atmospheric and Oceanic Science, University of Maryland, College Park, College Park, Marylandb
NOAA/Atlantic Oceanographic and Meteorological Laboratory, Miami, Florida

Monthly Weather Review (Proof Only)(Manuscript received 12 April 2022, in ﬁnal form 21 September 2022)
ABSTRACT: Particle ﬁlters avoid parametric estimates for Bayesian posterior densities, which alleviates Gaussian assumptions in nonlinear regimes. These methods, however, are more sensitive to sampling errors than Gaussian-based techniques such as ensemble Kalman ﬁlters. A recent study by the authors introduced an iterative strategy for particle ﬁltersthat match posterior moments}where iterations improve the ﬁlter’s ability to draw samples from non-Gaussian posterior
densities. The iterations follow from a factorization of particle weights, providing a natural framework for combining particle ﬁlters with alternative ﬁlters to mitigate the impact of sampling errors. The current study introduces a novel approach
to forming an adaptive hybrid data assimilation methodology, exploiting the theoretical strengths of nonparametric andparametric ﬁlters. At each data assimilation cycle, the iterative particle ﬁlter performs a sequence of updates while the priorsample distribution is non-Gaussian, then an ensemble Kalman ﬁlter provides the ﬁnal adjustment when Gaussian distributions for marginal quantities are detected. The method employs the Shapiro–Wilk test to determine when to make the transition between ﬁlter algorithms, which has outstanding power for detecting departures from normality. Experiments using
low-dimensional models demonstrate that the approach has a signiﬁcant value, especially for nonhomogeneous observation
networks and unknown model process errors. Moreover, hybrid factors are extended to consider marginals of more than
one collocated variables using a test for multivariate normality. Findings from this study motivate the use of the proposed
method for geophysical problems characterized by diverse observation networks and various dynamic instabilities, such as
numerical weather prediction models.
SIGNIFICANCE STATEMENT: Data assimilation statistically processes observation errors and model forecast
errors to provide optimal initial conditions for the forecast, playing a critical role in numerical weather forecasting. The
ensemble Kalman ﬁlter, which has been widely adopted and developed in many operational centers, assumes Gaussianityof the prior distribution and solves a linear system of equations, leading to bias in strong nonlinear regimes. On the other
hand, particle ﬁlters avoid many of those assumptions but are sensitive to sampling errors and are computationally expensive. We propose an adaptive hybrid strategy that combines their advantages and minimizes the disadvantages of
the two methods. The hybrid particle ﬁlter–ensemble Kalman ﬁlter is achieved with the Shapiro–Wilk test to detect theGaussianity of the ensemble members and determine the timing of the transition between these ﬁlter updates. Demonstrations in this study show that the proposed method is advantageous when observations are heterogeneous and whenthe model has an unknown bias. Furthermore, by extending the statistical hypothesis test to the test for multivariate
normality, we consider marginals of more than one collocated variable. These results encourage further testing for real
geophysical problems characterized by various dynamic instabilities, such as real numerical weather prediction models.
AU1

KEYWORDS: Kalman ﬁlters; Numerical analysis/modeling; Data assimilation; Nonlinear models

1. Introduction
For convection-permitting numerical weather prediction
systems, assimilating remotely sensed observation networks(e.g., radar and cloudy radiance measurements) is required to
depict mesoscale weather features accurately (e.g., Vukicevic
et al. 2004; Stengel et al. 2009; Privé et al. 2013). It is wellknown that strongly nonlinear model dynamics and observation operators, however, can induce bias in Gaussian-based
data assimilation methods that are commonly used for numerical weather prediction (e.g., Bocquet et al. 2010). Since ensemble Kalman ﬁlters (EnKFs; Evensen 1994; Houtekamer
and Mitchell 1998; Evensen and van Leeuwen 2000) approximate prior densities using a Gaussian and solve a linear system

Corresponding author: Kenta Kurosawa, kkurosaw@umd.edu

of equations to adjust a sample of model states to ﬁt the posterior mean and covariance, strongly nonlinear model dynamics
and measurement operators can lead to bias, which impedes
achieving accurate convection-permitting initial conditions fornext-generation weather forecast models. This limitation is apparent for multiscale weather prediction systems that exhibit
large uncertainty in smaller scales, or when observations are
sensitive to cloud processes (e.g., Poterjoy et al. 2017; Poterjoy
2022a). As a result, most infrared satellite data assimilation
studies mainly focus on clear-sky observations (e.g., Errico
et al. 2007; Fabry and Sun 2010; Minamide and Zhang 2017;
Honda et al. 2018). Therefore, developing new data assimilation methods that mitigate Gaussian assumptions is an active
area of research.
One strategy, which has gained momentum in recent years,
is to apply dimension-reduction procedures (viz., localization)
to particle ﬁlters (PFs; Penny and Miyoshi 2016; Poterjoy and

DOI: 10.1175/MWR-D-22-0108.1
Ó 2022 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright
Policy (www.ametsoc.org/PUBSReuseLicenses).

JOBNAME: MWR 00#00 2022 PAGE: 2 SESS: 8 OUTPUT: Sat Dec 24 08:35:27 2022 Total No. of Pages: 21/ams/mwr/0/mwrD220108

Monthly Weather Review (Proof Only)

2

MONTHLY WEATHER REVIEW

Anderson 2016; Poterjoy et al. 2017, 2019; Potthast et al.
2019). PFs avoid the parametric estimation of Bayesian posterior densities, thus providing great ﬂexibility for solving a
range of complex data assimilation problems. These methods,
however, are more sensitive to sampling errors than EnKFs.
As such, computational limitations pose a major obstacle,
which has limited research examining the potential of PFs for
operational weather prediction. Incorporating statistics froma large number of high-resolution ensemble members into thedata assimilation step is one of the most effective ways tomitigate the effects of sampling errors, but this strategy is notoften tractable.
Given the challenges discussed above, a natural progression
is to combine PFs with methods that rely on parametric density estimates when appropriate. Several papers have proposed to hybridize PFs with EnKFs (Stordal et al. 2011; Frei
and Kunsch 2013; Slivinski et al. 2015; Chustagulprom et al.
2016; Grooms and Robinson 2021) and with variational methods (Morzfeld et al. 2018). These methods are remarkably
accurate for cases of “moderate nonlinearity,” which are characteristic of situations with a non-Gaussian priors but Gaussianlike posterior distribution (Metref et al. 2014; Morzfeld andHodyss 2019; Grooms and Robinson 2021). For example, Frei
and Kunsch (2013) introduced a procedure that makes a continuous transition between the ensemble and the particle
ﬁlter update by factoring the likelihood. They choose a“splitting factor” to ensure that an effective ensemble size
is maintained within a certain tolerance of a user-speciﬁed
threshold. Based on this approach, Chustagulprom et al.
(2016) developed a method to hybridize the general linear
ensemble transform ﬁlter (LETF) framework and PFs, which
can use observation-space localization and avoid linear assumptions for observation operators. Grooms and Robinson(2021) also introduced a ﬁlter that combines PFs with EnKFs,
which is generally similar to Frei and Kunsch (2013) andChustagulprom et al. (2016) in that it factors the likelihood.
This method, like others, is effective for problems characterized by medium nonlinearity in model dynamics or measurement operators. In these papers, the value of the splitting
factor is determined adaptively by the effective ensemble size.
This choice is a heuristic one, thus motivating additional research into how to optimally combine PFs with EnKFs. For
example, Nerger (2022) propose a method for estimating hybrid coefﬁcients that is based not only the effective ensemble
size but also the kurtosis and skewness of the ensemble. Thismethod still requires the tuning of hyper-parameters, but
generates more accurate ﬁlter estimates than using effective
ensemble size alone. As mentioned in Chustagulprom et al.
(2016), a more powerful and computationally feasible alternative is to adopt Kullback–Leibler divergence (KL divergence; Kullback and Leibler 1951) as a means of identifying
proper choices of prior error distribution. This approach isone of the most frequently used objective functions to measure deviations from Gaussianity in forecast error distribution for weather models (Kondo and Miyoshi 2019; Li et al.
2019; Ruiz et al. 2021; Pimentel and Qranfal 2021). However, it is difﬁcult to measure non-Gaussianity by the KL divergence for high-dimensional systems when the ensemble

VOLUME 00

size is small or when a strange attractor makes numericalconvergence and proper deﬁnition of the continuous limit
complicated (Bocquet et al. 2010).
In this study, we introduce a novel approach to forming an
adaptive PF–EnKF data assimilation methodology, which exploits the theoretical strengths of nonparametric (PF) andparametric (EnKF) ﬁlters. For this purpose, we use a recently
proposed PF (Poterjoy 2022a,b), which introduces an iterative
strategy for PFs that match posterior moments. For this PF,
iterations improve the ﬁlter’s ability to draw samples from
non-Gaussian posterior densities despite ﬁtting a limited number of moments. The iterations follow from a factorization ofparticle weights, which also provide a straightforward meansof combining PFs with EnKFs to reduce the impact of sampling errors. To achieve the adaptive mixed methodology ateach data assimilation cycle, we repeat the iterative PF update
while the prior sample distribution is non-Gaussian, and update with EnKF when Gaussian distributions for marginal
quantities are detected. Here, we introduce a statistical hypothesis testing approach to determine when to make thistransition between ﬁlter algorithms. Several papers on data assimilation use statistical hypothesis tests for normality to measure the difference between the prior or posterior distribution
and the normal distribution (e.g., Bocquet et al. 2010; Poterjoy2016). However, research of incorporating a normality test directly into assimilation processes is unexplored.
The current study ﬁrst compares the power of several hypothesis tests by performing Monte Carlo simulations of data
generated from choice distributions that are often used tocharacterize errors. Then we examine the newly developed
mixed hybrid methodology employing the Shapiro–Wilk test
(Shapiro and Wilk 1965), which has outstanding power among
omnibus tests for detecting departures from normality (e.g.,
Srivastava and Hui 1987; Mendes and Pala 2003; Farrell et al.
2007; Villaseñor and González-Estrada 2009). Using the hypothesis testing approach allows us to accurately detect Gaussianity, even with small ensemble sizes, and explore possibilities
other than effective ensemble size and the KL divergence fordetermining the splitting factor adaptively.
The manuscript is organized in the following manner. In
section 2, we brieﬂy review the four well-known normality
tests and compare the power of the tests. Section 3 introduces
a statistical hypothesis testing approach to forming an adaptive PF–EnKF hybrid. We discuss the results and ﬁndings ofnumerical experiments conducted using low-dimensional toy
models in section 4. The last section discusses major ﬁndings
from this study.

2. Power comparisons of four normality tests
The assumption of a normal distribution is often an underlying premise of many academic ﬁelds and studies, including
data assimilation. When the assumption of normality is violated, interpretations and inferences may lack reliability andvalidity. There are three commonly used procedures for evaluating whether a random independent sample comes from anormal population: graphical methods (histograms, box plots,
Q–Q plots), moment-based methods (skewness and kurtosis),

JOBNAME: MWR 00#00 2022 PAGE: 3 SESS: 8 OUTPUT: Sat Dec 24 08:35:27 2022 Total No. of Pages: 21/ams/mwr/0/mwrD220108

MONTH 2022

and formal normality tests. In particular, a signiﬁcant amountof normality tests have been proposed, and several studies
have already compared their power (e.g., Dufour et al. 1998;
Thadewald and Büning 2007; Razali and Wah 2011; Saculinggan
and Balase 2013). In this section, we compare the power offour well-known formal tests of normality: Shapiro–Wilk test
(SWT; Shapiro and Wilk 1965), Kolmogorov–Smirnov test
(Kolmogorov 1933; Smirnov 1939), Lilliefors test (Lilliefors
1967), and Anderson–Darling test (Anderson and Darling1954). The following subsections brieﬂy review the four normality tests, and describe the simulation procedure.

a. Methodology for the four hypothesis tests
1) SHAPIRO–WILK TESTMonthly Weather Review (Proof Only)

3

KUROSAWA AND POTERJOY

The normality test introduced by Shapiro and Wilk (1965)
is the ﬁrst test to detect departures from normality with skewness, kurtosis, or both (Althouse et al. 1998). It has become
the most potent omnibus test in most situations because of its
good power properties compared to a wide range of alternative
tests. The basic idea behind SWT is to measure the goodness of
ﬁt of a straight line to a normal Q–Q plot (linear regression).
Given an ordered random sample, x1 , x2 , … , xn, the original SWT statistic is deﬁned as!
"2n∑ ai xi

W5 n

i51

,

Kolmogorov–Smirnov test was ﬁrst proposed by Kolmogorov(1933) and then improved by Smirnov (1939). The one-sample
Kolmogorov–Smirnov test is a nonparametric test of the null
hypothesis that the population cumulative distribution function(cdf) of the data is equal to the hypothesized cdf. Given an ordered random sample, the statistic is deﬁned as
D 5 max[|F * (x) 2 Fn (x)|],x3) LILLIEFORS TEST
Kolmogorov–Smirnov test is appropriate when the hypothesized distribution parameters are completely known because
the null distribution must be completely speciﬁed. In contrast,
Lilliefors test, which is a modiﬁcation of Kolmogorov–Smirnov
test introduced by Lilliefors (1967), is a goodness-of-ﬁt test forsituations where the parameters of the null distribution are unknown and have to be estimated. Given an ordered random
sample, the Lilliefors test statistic is deﬁned as
D 5 max[|F * (x) 2 Sn (x)|],
x

2

i51

where x is the sample mean, and ai is the expected values of theorder statistics of independent and identically distributed random variables sampled from the standard normal distribution:
mT V21(mT V21 V21 m)1/2,(2)

where m 5 (m1 , …, mn )T . Here, the vector m consists of the
expected values of the order statistics of independent randomvariables with identical distributions that are sampled from anormal distribution, and V is the covariance matrix of thoseorder statistics.
The null hypothesis of SWT is that the data originate froma normally distributed population. Small values of W lead tothe rejection of the null hypothesis, where 0 # W # 1. The
original SWT was limited to a sample size of 50 or less. Royston(1982) extended SWT to large samples and provided an approximation of the test statistic W, and Royston (1983) suggested a
test for multivariate normality based on SWT. Royston (1992)
arrived at an improved approximation to the weights, which allows SWT to effectively detect departures from multivariate
normality for smaller sample sizes. Last, Royston (1995) introduced the FORTRAN algorithm AS R94, which is used in thecurrent study. The algorithm includes a scaling process that sets
the mean of the sample to zero and a centering process that sets
the variance of the sample to one. Therefore, the algorithm allows the results to generalize not only to the standard normal
distribution, but also to a normal distribution with a nonzero
mean and nonunit variance.

(3)

where F*(x) is the cdf of the hypothesized distribution, andFn(x) is the empirical cdf. When the statistic value D is signiﬁcant, the hypothesis that the sample comes from a normally
distributed population is rejected.

(1)∑ (xi 2 x)

ai 5 (a1 , …, an ) 5

2) KOLMOGOROV–SMIRNOV TEST

(4)

where F*(x) is the cdf of the hypothesized distribution, andSn(x) is the empirical cdf. The Lilliefors test statistic is thesame as the Kolmogorov–Smirnov test statistic, but the tablesof critical values of the two tests are different, leading to different conclusions and decisions. In Kolmogorov–Smirnov
test, we must completely give the null distribution. On theother hand, Lilliefors test is a two-sided goodness-of-ﬁt test
and is powerful when the parameters of the null distribution
are unknown.

4) ANDERSON–DARLING TEST
Anderson–Darling test, introduced by Anderson andDarling (1954), is a modiﬁcation of the Cramer–von Miles
test (Cramér 1928). The Anderson–Darling test statistic is
deﬁned as#‘n[Fn (x) 2 F(x)]2 w(x)dF(x),(5)
2‘where n is the sample size, w(x) is a weight function, F(x) is
the hypothesized distribution, and Fn(x) is the empirical cdf.
The weight function is deﬁned asw(x) 5 {F(x)[1 2 F(x)]}21 ,(6)

and Arshad et al. (2003) suggested the following formula as
the test statistic of Anderson–Darling test:n2i 2 1{ln[F(Xi )] 1 ln[1 2 F(Xn112i )]},ni51

A2n 5 2 n 2 ∑(7)

JOBNAME: MWR 00#00 2022 PAGE: 4 SESS: 8 OUTPUT: Sat Dec 24 08:35:28 2022 Total No. of Pages: 21/ams/mwr/0/mwrD220108

4

MONTHLY WEATHER REVIEWMonthly Weather Review (Proof Only)

TABLE 1. Classiﬁcation of cases by skewness and kurtosis of the
distribution.
%$$$$
Skewness ( b1 ),kurtosis (b2)
CaseDistributions%$$$$
b1 5 0, b2 5 3N(0, 1)A%$$$$
b 5 0, b2 , 3U(0, 1), Beta(2, 2)
B%$$$1$
b1 5 0, b2 . 3
Logistic(0, 1), t (5)C%$$$$
b1 Þ 0, b2 , 3
Weibull(3, 5), Beta(2, 4)
D
%$$$$
E
b1 Þ 0, b2 . 3
Gamma(2, 2), x2(10)where X1 , X2 , … , Xn are the ordered sample data
points. The weighting function (6) is more sensitive to outliers
because the weights of the observed values at the edges of the
distribution are larger. Therefore, it is especially suited for detecting deviations from normality at the tails of the distribution.

b. Simulation procedures

T1

Monte Carlo simulations are the most commonly used approaches to compare and evaluate the accuracy of a hypothesis test in detecting the degree of contamination by outliers
and the test power with respect to sample size. Following previous studies, we use Monte Carlo simulations to evaluate thepower of SWT, Kolmogorov–Smirnov test, Lilliefors test, andAnderson–Darling test statistics in testing whether a random
sample of n independent observations is obtained from a population with a normal N(m, s2) distribution. Several papers
have already shown the superiority of SWT over the othertests (e.g., Mendes and Pala 2003; Razali and Wah 2011).
Thus, the simulation in this section focuses on a larger selection of distributions than previous studies, which are motivated by the diverse shapes of error distributions found forgeophysical models and observing systems, and aims to reconﬁrm the superiority of SWT over the other tests.
As summarized in Table 1, we examine the following ﬁve
cases of nine distributions to cover a variety of standardized%$$$$
skewness ( b1 ) and kurtosis (b2): N(0, 1), U(0, 1), Beta(2, 2),Logistic(0, 1), t(5), Weibull(3, 5), Beta(2, 4), Gamma(0, 1),%$$$$
and g2(10). The values of b1 and b2 for each distribution aresummarized in Table 1. For each distribution, we set the signiﬁcance level at 0.05, the sample sizes at n 5 5, 20, 40, 100,
300, 500, 1000, and the number of trials at 100 000. The nulland alternative hypotheses of the four tests are as follows:
H0 : The distribution is normal
H1 : The distribution is not normal:

T2

As summarized in Table 2, the “test power” (true positive) ofa hypothesis test is the probability that the test correctly rejects the null hypothesis when the alternative hypothesis istrue. On the other hand, the “type I error” (false positive) isthe error of rejecting the null hypothesis when it is actually true.
Therefore, if a sample is taken from N(0, 1) population, thenumber of rejected H0 hypotheses is the probability of a type I
error (case A in Table 1). In contrast, if the samples are from a

VOLUME 00

TABLE 2. Deﬁnitions of terminologies in a statistical test.
Decision

Null hypothesis
H0 is

True

False

Accept H0

Reject H0

Speciﬁcity“True negative”
Probability: 1 2 a
Type II error“False negative”
Probability: b

Type I error“False positive”
Probability: aTest power“True positive”
Probability: 1 2 b

population that is not normal distribution, the number of H0 rejected is the power of the test (cases B–E in Table 1).
In Table 2, the probability of a type I error occurring isdenoted by a and the probability of a type II error by b. The
probability of both the type I and type II errors should besmall. However, it is impossible to make both small because
the risk ratios a and b are in a trade-off relationship. Generally speaking (in society), committing a type I error is often
a more serious problem. Moreover, as discussed in detail in
section 3, in the current study, we prefer erring on the side ofa lower type I error in order to avoid the use of an EnKF
when the distribution is truly non-Gaussian. This is because
the PF update can provide adequate estimates for both nonGaussian or Gaussian prior distributions; the same cannot be
said about the EnKF. Therefore, the correct procedure for hypothesis testing is to determine the acceptable risk rate a inadvance, and then select the hypothesis test method with thehighest test power 1 2 b among them. Hence, in the current
study, we focus on the type I error for N(0, 1) and on the test
power for the other distributions (cases B–E in Table 1); We
will not look at the type II error and speciﬁcity in this study.

c. Results
Figure 1 shows the variation of the type I error (Fig. 1a)and test power (Figs. 1b–i) with the sample sizes n for thefour tests for each distribution when a 5 0.05. When the distribution is N(0, 1), all four tests generally achieved a 5 0.05.
%$$$$
In the case of the symmetric distributions ( b1 5 0 ; Figs. 1b–e),
SWT is the best, followed by Anderson–Darling test, Lilliefors
test, and Kolmogorov–Smirnov test. However, all tests have low
power when the sample size is less than 100. In particular, whenb2 is greater than 3, the power of Kolmogorov–Smirnov test issigniﬁcantly inferior to the other three (Figs. 1d,e). All othertests can attain 80% power when the sample size is 1000. The%$$$$
power for the asymmetric distributions ( b1 Þ 0 ; Figs. 1f–i) is
also highest for SWT, followed by Anderson–Darling test, Lilliefors test, and Kolmogorov–Smirnov test. For the Weibull(3, 5)
distribution, the overall power is low because the kurtosis b2 isclose to 3 (Fig. 1f). In the other cases, however, (Figs. 1g–i),
SWT and Anderson–Darling test require 200 samples to achieve
90% power, while Kolmogorov–Smirnov test requires 500 samples; Lilliefors test shows power for sample sizes between these
two numbers.
The overall results show that in all cases the power of SWT
is superior to the other tests for small sample sizes, which is

F1

JOBNAME: MWR 00#00 2022 PAGE: 5 SESS: 8 OUTPUT: Sat Dec 24 08:35:28 2022 Total No. of Pages: 21/ams/mwr/0/mwrD220108

Monthly Weather Review (Proof Only)

MONTH 2022

KUROSAWA AND POTERJOY

5

FIG. 1. (a) Type I error and (b)–(i) test power of Shapiro–Wilk test (SW; orange), Anderson–Darling test (AD; cyan), Lilliefors test
(LF; blue), and Kolmogorov–Smirnov test (KS; green) for different distributions and sample sizes. The magenta line shows each distribution, and thedashed line shows a close normal distribution to each distribution. The distributions cover a variety of standardized%$$$black$skewness ( b1 ) and kurtosis (b2).

the regime of interest for ensemble data assimilation applied toweather models. This result is generally consistent with theresults of many other previous studies mentioned before.
Therefore, this paper uses SWT to detect Gaussianity in the
new hybrid method hereafter. Note that “failing to reject thenull hypothesis” is not the same as “accepting the null hypothesis.” In such a case, it is still not exactly clear whether the nullor alternative hypothesis is correct. For simplicity, this study
interprets the null hypothesis to be that “the samples are froma population that follows a normal distribution.” If the null hypothesis of SWT is not rejected then the samples are assumed
to come from a Gaussian distribution.

3. Implementation with the local particle filter
This section presents the mathematical framework for implementing the adaptive PF–EnKF hybrid method by embedding SWT, the most powerful statistical test in the previous
section, into the local PF. For the purposes of the adaptive
hybrid methodology, the current study uses the recently
proposed PF by Poterjoy (2022b, hereafter P22), which introduces an iterative strategy. We brieﬂy describe the partsof PF that are relevant to the implementation of the proposed method.
The local PF operates by assimilating observations withindependent errors sequentially and combining sampled

Fig(s). 1 live 4/C

particles and prior particles for each observation. By serially processing an observation y in a sequence of observations and updating particles after each observation space
sampling step, posterior particles can be adjusted in a mannerconsistent with bootstrap sampling. The nth updated particle
xny is represented by the linear combination of the resampled
particle xkn , conditioned on all observations before y, and theprior particle xn as follows:
xny 5 x y 1 r1 + (xkn 2 x y ) 1 r2 + (xn 2 x y ),(8)

where kn is the index of each sampled particle, and x y is the
localized posterior mean accumulating the full weight of allobservations up to y. r1 and r2 are vectors of coefﬁcients that
ensure the update satisﬁes the posterior mean and variance ofmarginals everywhere in state space}as depicted by importance weights.
Poterjoy et al. (2019) introduced several ﬁlter stabilization
strategies in the PF of Poterjoy (2016) to avoid particle degeneracy. In particular, regularization and tempering are effective
methods when sampling errors are large, and the sample size
is small. Regularization is equivalent to increasing the particle
weights up to a power b by inﬂating the observation error
variance. This regularization allows the particles to acquirea speciﬁc “effective sample size Neff,” and is particularly
helpful in stabilizing the ﬁlter when all particles are far

JOBNAME: MWR 00#00 2022 PAGE: 6 SESS: 8 OUTPUT: Sat Dec 24 08:35:28 2022 Total No. of Pages: 21/ams/mwr/0/mwrD220108

Monthly Weather Review (Proof Only)

6

MONTHLY WEATHER REVIEW

from an observation. The regularization provides temporary iterations for the local PF, which is a posterior tempering strategy. This iterative approach also improves the ﬁlter’sability to sample from non-Gaussian posterior densities,
even though it ﬁts a limited number of moments. The iterations consist of a factorization of particle weights, thus providing a natural framework for combining the local PF withalternative ﬁlters to reduce the impact of sampling error.
When provided with Gaussian likelihoods, a partial update
performed by a PF can adjust particles to more closely resemble samples from a Gaussian, even if the prior exhibits
a complex nn-Gaussian shape. The resulting intermediate
update then makes the EnKF an appropriate choice for theremaining update (Grooms and Robinson 2021). P22 introduced the hybrid parameter k, which is an Nx-dimensional
vector that determines when to switch from a PF update toa parametric ﬁlter update. For the iterative PF, the hybridt” areparameter k and “the target effective sample size Neff
required to be speciﬁed by users. PF updates are repeated untilNk∑k51 bj,k 5 kj for 0 # kj # 1 at the jth grid point, where Nk

is the number of iterations. Here, bj enforces a minimum
Neff for weight at the jth grid point. When Neff is belowt
Neff, b is determined adaptively by Eq. (28) in Poterjoyt
et al. (2019) so that Neffis satisﬁed. Following the initial setof local PF iterations, the last adjustment is performed using an EnKF with the measurement error variance R inﬂated by the factor 1/hj , where hj 5 1 2 kj. For example, tohybridize the PF and EnKF in the ratio of 7:3 at the jth grid
point, kj is set to 0.7 in the ﬁrst place, and PF updates areNkrepeated until ∑k51
bj,k 5 kj 5 0:7. The value of bj,k is det
in the kth iteration,
termined adaptively based on Neffwhich means Nk is determined adaptively as well and different at each grid point. Therefore, since Nk at jth gridtispoint is determined by bj,k, Nk becomes larger when Neff
set to a larger value, and vice versa. Last, an EnKF update
is performed with R inﬂated by a remaining factor hj, where
hj 5 1 2 kj 5 0.3. Note that, for simplicity, P22 uses thesame value for kj at all grid points and in all data assimilation cycles.
In the current study, we allow k and h to be adjusted
adaptively through space and time during data assimilation,
while these values are held constant and set heuristicallythrough tuning in P22. We repeat the PF update until SWT
Nk
bj,k upsuggests that particles are samples from a Gaussian; ∑k51
to this point is deﬁned as kj, and the value of hj is determinedby 1 2 kj. Once kj at all grid points has been determined, we
perform the serial ensemble square root ﬁlter (serial EnSRF;
Whitaker and Hamill 2002) with R inﬂated by the inverse of
hj as the ﬁnal adjustment. Thus, the Kalman gain matrix at
the jth grid point when the ith observation is assimilated is
described as follows:
!
"21
1
f fT
f fT
(9)
Kj 5 Ej Di Di Di 1 Ri ,hjwhere Ef consists of model-space forecast ensemble perturbations and Df consists of observation-space forecast ensemble

VOLUME 00√$$$$$$$$$$
perturbations, with both matrices normalized by 1/ Ne 2 1.
Note that the use of a tangent linear measurement operatorin (9) is avoided in the current study and most others by assimilating observations serially.
In the case that SWT does not detect Gaussianity during Nk
iterations at the jth grid point, kj and hj become 1 and 0, respectively, and no EnKF update is performed at the grid
point. Similarly, if SWT detects Gaussianity in the ﬁrst iteration at the jth grid point, then kj 5 0 and hj 5 1, and no PF
update is performed at the grid point. Thus, in situations
where the posterior is clearly non-Gaussian, the ﬁlter canhave the option of retaining the local PF. The hybrid approach aims to obtain an intermediate distribution that iscloser to Gaussian than the prior distribution by the PF updates, and to make this intermediate distribution closer toGaussian by the EnKF. In cases that we cannot obtain an intermediate distribution closer to Gaussian, we can performthe iterative PF updates alone, without using EnKF in the laststep, which is the strength of the adaptive strategy. The advantage of using the EnKF in the last step if a Gaussian is encountered during iterations is purely due to it being a morerobust choice when ensemble sizes are small (and the distribution is indeed Gaussian).
tstill needs to be speciﬁed by the user
As in P22, the Nefft
and this parameter can inﬂuence the results. In general, Neffdetermines when ﬁlter updates are made during iterations.
ttypically leads to more iterations and a larger ﬁnal
High Neffteffective ensemble size than a small Neff
. This choice is ultimately a trade-off between the frequency of performing SWT
and cost of implementation.
Since the computational cost of SWT is not expensive, theadaptive approach introduced in the current study, which incorporates the statistical test into the local PF introduced by
P22, is generally less computationally expensive than the iterative LPF. Under most circumstances, the hybrid requires
fewer iterations, thus leading to a cost saving. Nevertheless,the PF introduced by P22 is computationally more costly thanpure EnKF because of the use of regularization and tempering. For more information, please refer to Poterjoy (2022a)
and P22.
Note that k and h are uniquely speciﬁed for each observationspace prior variable as well. In this case, they are Nydimensional vectors and we again use SWT to determine
when each element of observation-space forecast ensembles
may follow a Gaussian distribution. The k and h deﬁned forthe observation-space are used for the observation-spaceﬁlter updates.
In summary, the adaptive hybrid PF–Serial EnSRF with
SWT are realized by Algorithms 1–2. In both algorithms, xf
and xa are Nx-dimensional background and analysis vectors,
respectively, and yo is an Ny-dimensional set of observations.
The observation operator H maps a model state to its corresponding observation state:
yf 5 Hi (xf ),(10)

where Hi is the measurement operator for the ith observation.

JOBNAME: MWR 00#00 2022 PAGE: 7 SESS: 8 OUTPUT: Sat Dec 24 08:35:29 2022 Total No. of Pages: 21/ams/mwr/0/mwrD220108

MONTH 2022
AU2

7

KUROSAWA AND POTERJOY

Algorithm 1 Adaptive mixed PF–EnKF update with SWT
1: function pf_enkf_hybrid
2:

k51

3:

k 5 0.04:k_residual 5 1.0 – k

5:while max(k_residual) . 0 do

. Vector with Nx dimensions

if k_residual (j) . 0 then7:f

SWT_result $ WT(xj )8:. Shapiro-Wilk Test

if SWT_result 5 Gaussian then9:10:

hj 5 k_residual(j)k_residual(j) 5 0.0

11:Monthly Weather Review (Proof Only). Tempering

for j 5 1: Nx do

6:

end if12:
end if

13:
14:

end for15:(bk, k_residual) $ Regularization(k_residual)

16:

for i 5 1: Ny do

17:

xa $ The Local PF (xf , yoi , bk )18:

xf $ xa

19:

end for

20:

k5k11

21:

end while22:

xa $ EnKF_tempered(xf, yo, h). The Local PF core

. EnKF as the last adjustment

23: end function

4. Numerical experiments with low-order models

Algorithm 2 Serial EnSRF update with inﬂated Rfo

1: function EnKF_tempered(x , y , h)2:

for I 5 1: Ny do
yf 5 Hi (xf )

3:

1
f
f
Ef 5 %$$$$$$$$$$ [d x1 | · · · |d xN ]e
Ne 2 14:

1f
f
f
Di 5 %$$$$$$$$$$ [d y1 | · · · |d yN ]e
Ne 2 1

5:for j 5 1: Nx do

6:

!

1
f fT
f fT
Kj 5 Ej Di Di Di 1 Ri
hj

7:
8:

end for9:11:

x a 5 x f 1 K(yoi 2 yf )&$$$$$$$$$$$$$$$$21
Ri


a 5 1 1
f fT
Di Di 1 Ri

12:

Ea 5 Ef 2 K̃Di10:

K̃ 5 aK

f

a

a

a

13:

x 5x 1E14:

xf $ xa

15:

end for

16: end function"21

This section explores the behavior of the newly developed
method through numerical simulations. In the ﬁrst experiment, we use a simple univariate problem to illustrate the difference of the adaptive hybrid method between iterative
EnKF and bootstrap PF using tempering. In the second experiment, we use the 40-variable dynamic model of Lorenz(1996) to compare the advantages of the adaptive method
over EnKF, the local PF, and hybrid PF–EnKF with ﬁxed values of k and h. These experiments use simulated measurements to target several scenarios, such as varying spatial
density, highly nonlinear dynamics, mixed measurement operators, and unresolved model error. The last experiment uses
an idealized kinematic vortex, which was used in Poterjoy(2022a) to replicate ﬁndings from real-data applications. The
kinematic vortex model allows us to generate observations
that emulate realistic observations for an application containing large spatial dependence across variables, while
retaining great ﬂexibility in our construction of data assimilation experiments. Among the common parameters used
for idealized data assimilation applications, such as observation error variance, observation locations, and ensemble size,these experiments contain parameters that indirectly control
the shape of the full multivariate prior, thus allowing for ananalysis of the adaptive hybrid technique under controllable
conditions.

JOBNAME: MWR 00#00 2022 PAGE: 8 SESS: 8 OUTPUT: Sat Dec 24 08:35:29 2022 Total No. of Pages: 21/ams/mwr/0/mwrD220108

Monthly Weather Review (Proof Only)

8

MONTHLY WEATHER REVIEWVOLUME 00

FIG. 2. A univariate example of how the updates differ in each iterative ﬁlter for (left ) EnKF–MDA, (center) IPF, and (right) adaptive
mixed IPF–EnKF. Each row corresponds to an iteration. The blue and red lines indicate marginal prior and posterior pdfs, respectively.
The black dashed line indicates the observation likelihood.

a. Univariate application
Using a univariate example, we can visualize how the newly
proposed adaptive hybrid method works compared to ﬁltersthat use iterative strategies. We compare three iterative ﬁltersin this section: EnKF with the multiple data assimilation
scheme (EnKF–MDA) proposed by Emerick and Reynolds(2012), bootstrap ﬁlter adopting the iterative approach (IPF),
and a hybrid of the IPF and EnKF (adaptive IPF–EnKF). The
number of iteration is set to four for EnKF–MDA and IPF
for this demonstration. For EnKF–MDA, when the same observation is assimilated Na times, the inﬂated measurement
error covariance matrix is used:
K 5 Ef Df T (Df Df T 1 ai R)21 ,

(11)

where
Na

∑

1

i51 ai

5 1:(12)

In this experiment, we use ai 5 4 for i 5 1, … , Na, where
Na 5 4. For further details on EnKF–MDA, we encourage
readers to review the mathematical descriptions in Emerick andReynolds (2012). The IPF also uses a factorization of the likelihood to break the PF update step into a sequence of fourFig(s). 2 live 4/C

updates, namely, bi 5 1/4 for i 5 1, … , Na, where Na 5 4. For
the adaptive IPF–EnKF, we ﬁrst set bi 5 1/8 and repeat thebootstrap PF update until SWT detects that prior members are
samples from a Gaussian. We then replace the remaining PF
update with an EnKF update using R inﬂated with the inverseof the remaining likelihood [1/(1 2 k) 5 1/h]. Note that themultiple updates in EnKF–MDA and IPF are identical to single
updates of each because the operator is linear in this example.
Consider the example shown in Fig. 2, where 104 priormembers are updated using an observation whose value is 4
and observation error standard deviation is set to sy 5 0.8.
Among the 104 prior members, three quarters are selectedfrom N(24, 1.22), while the rest are from N(3.5, 1.22), whosemean value is smaller than the observation. Therefore, theprior ensemble is a bimodal distribution. Using the same prior
for all three ﬁlters, Figs. 2a–c show the posterior distribution
after the ﬁrst iteration. In all cases, we can see that each ﬁlter
shifts the ensemble toward the observation; however, the
EnKF–MDA inherits the bimodal distribution of the prior forthe posterior distribution, while the IPF and IPF–EnKF correctly
retain a single mode.1 The bimodal posterior distributions in the1
For the provided ensemble size, the last iteration of the IPF isan accurate estimate of the true Bayesian posterior.

F2

JOBNAME: MWR 00#00 2022 PAGE: 9 SESS: 8 OUTPUT: Sat Dec 24 08:35:29 2022 Total No. of Pages: 21/ams/mwr/0/mwrD220108

9

KUROSAWA AND POTERJOYMONTH 2022

TABLE 3. Conﬁguration of cycling data assimilation experiments.
Experimental settings
Expt name

H(x)F

Observation points

Linear Case
Nonlinear Case 1
Nonlinear Case 2
Mix Case 1
Mix Case 2
Mix Case 3xx2
log(|x|)x and x2
x and log(|x|)
x2 and log(|x|)

8.0; 9.0
8.0; 9.0
8.0; 9.0
8.0
8.0
8.0

Evenly distributed; data void
Evenly distributed; data void
Evenly distributed; data void
Evenly distributed
Evenly distributed
Evenly distributed

Monthly Weather Review (Proof Only)

Parameter filter settings
Ne

ROI

a, g

10
20
40
100
300

2; 1
5; 1
7; 2
9; 3
9; 3

0.3
0.3
0.3
0.3
0.3

EnKF–MDA are not relieved by the completion of all iterations(Fig. 2j). In the IPF, after all the updates, the posterior pdf is relatively close to the likelihood of the observation, but exhibits
negative skewness because numerous particles remain in theleftmost mode (Fig. 2k). In the adaptive IPF–EnKF case,
SWT detected Gaussianity in the distribution of the ensembles after three iterations of the bootstrap PF (Fig. 2i), andthen EnKF was performed using R inﬂated by 8/5, which is
the inverse of the remaining observation error variance (Fig. 2l).
As a result, the IPF–EnKF posterior is close to the IPF, indicating that the hybrid method correctly transitioned to the partial
EnKF step once a Gaussian distribution was detected. Furthermore, we emphasize that the univariate application is presented
for illustration only, as the IPF–EnKF is not expected to provide
beneﬁts over the IPF when the ensemble size is large.

b. 40-variable dynamical system
1) EXPERIMENTAL DESIGNSFor the next set of experiments, we assess the proposed
adaptive hybrid strategy through idealized numerical experiments with the Lorenz 40-variable model (Lorenz 1996;
Lorenz and Emanuel 1998), denoted L96 hereafter. The
model consists of Nx equally spaced variables and is deﬁned bydxi
5 (xi11 2 xi22 )xi21 2 xi 1 F,
dt

(13)

where i 5 1, 2, …, Nx with cyclic boundaries: xi1Nx 5 xi andxi2Nx 5 xi. The model is integrated forward numerically using
the fourth-order Runge–Kutta scheme and a model time stepof 0.05 nondimensional units, which is corresponding to 6 h(Lorenz 1996). As in Lorenz (1996), we ﬁx Nx at 40 and use
F 5 8.0, except for one set of experiments that consider an imperfect model case; in this case, measurements are simulated
from a model trajectory with F 5 8.0, but the model forcing F
is ﬁxed at 9.0.
Experiments include three forms of measurement operator
H as in Kurosawa and Poterjoy (2021): “Linear Case,”^,“Nonlinear Case 1,” and “Nonlinear Case 2” use H(x) 5 x^+x^ , and H(x) 5 log[ABS(^^
H(x) 5 xx )], respectively. Here, x
is a subset of Ny variables in x chosen by H, and ABS stands
for the absolute value of each element. Uncorrelated Gaussian
errors selected from N(0, s2y I) are added to each operator:
sy 5 1.0 for the ﬁrst two experiments, while sy 5 0.1 for thethird case because of the smaller information content provided by this observation. All experiments use Ny 5 20 observations applying one or two of the operators. When only one
observation operator is used, there are three settings: a settingwith evenly distributed observations and F 5 8.0 (“normal”),
a setting with evenly distributed but F 5 9.0 (“model error”),
and a setting with missing observations in some places (“datavoid”). In contrast to the setting where observations are homogeneous throughout the domain, the “data void” setting
is designed to target the heterogeneous observation network of real atmospheric models. Note that for this setting,
we set the observation points at grid points 1–10, 21–30. In
the “mix” case, two observation operators are used, namely,
the ﬁrst and the second half of the observation points use
different observation operators. The experimental settingsfor each of these cases are summarized in the section of
Table 3.
All experiments in this section use an observation frequency of 6 h. Observations are assimilated over a 10-yr period, and root-mean-square errors (RMSEs) from the last
9 years are used to quantify the accuracy of the posterior analyses, ignoring the ﬁrst year spin up period. In this set of experiments, we perform 100 parallel trials out of an abundanceof caution. For localization, we use the ﬁfth-order correlation
function controlled by a radius of inﬂuence (ROI) given by
Gaspari and Cohn (1999). For posterior inﬂation, the current
study adopts the strategy known as relaxation to prior perturbation (RTPP; Zhang et al. 2004) after EnKF update. Similar
to the a used in the relaxation method, for the local PF, we
use a mixing parameter g to maintain particle diversity during
updates in (8). When the ensemble size is small, this parameter works to prevent ﬁlter divergence. g is a scalar between

T3

JOBNAME: MWR 00#00 2022 PAGE: 10 SESS: 8 OUTPUT: Sat Dec 24 08:35:30 2022 Total No. of Pages: 21/ams/mwr/0/mwrD220108

Monthly Weather Review (Proof Only)

10

MONTHLY WEATHER REVIEWVOLUME 00

FIG. 3. Mean analysis RMSEs of the 12 experiments with different settings of k as a functionof ensemble size. Results are shown for (a) Linear Case, (b) Nonlinear Case1, and (c) Nonlinear
Case2. Values are from the average of the last 9 years with 100 parallel trials.

0 and 1, and acts to increase diversity in particles without modifying prior or posterior error variance. Each time the particles
in state space are updated, the prior particles are mixed withthe resampled particles (P22). The target Neff is ﬁxed att
5 0:5Ne for all experiments. We conﬁgure the number ofNeff
members 10, 20, 40, 100, and 300, and arbitrarily tune all ﬁlter
parameters, namely, ROI, a, and g, for each ensemble size.
For the “normal” and “mix” settings, ROIs for Ne 5 10, 20,
40, 100, and 300 are 2, 5, 7, 9, and 9, respectively. On the other
hand, to ensure the stability of the experiments, for the “datavoid” and “model error” settings, ROIs for Ne 5 10, 20, 40,
100, and 300 are 1, 1, 2, 3, and 3, respectively. The settings ofﬁlter parameters are summarized in the bottom section of
Table 3. Under each experimental setting, we performed a total of 12 experiments, one in which the value of k is estimatedadaptively, and the others in which k is ﬁxed at 0.1 incrementsfrom 0 to 1. Note that, in this section, the inﬂation and localization parameters for EnKF, the local PF, and hybrid experiments are uniﬁed, so we limited this tuning to experimentsthat use the LPF and EnKF alone. The tuning step is complicated for hybrid implementations, since we would have different optimal values for ROI and other parameters as soon as
we change the k value. This feature makes it difﬁcult to identify optimal parameters in a cost-effective manner. While we
acknowledge this limitation in the comparisons, we note thathybrid conﬁgurations still tend to outperform the LPF and
EnKF despite not following a rigorous tuning. In other words,
we believe that the use of optimal parameters may slightly
change the results of the following experiments, but it will notchange the conclusion of this section.

Fig(s). 3 live 4/C

2) RESULTS
We summarize the results for the “normal” setting in Fig. 3,
and the mean values of k in the adaptive experiments withthis setting are shown later (see Fig. 8a). In Linear Case
(Fig. 3a), when the number of members is small, the higher
the ratio of the EnKF, or the closer the value of k is to zero,
the lower the RMSE. However, with 40 members, the experiment that performs partial local PF is optimal (k 5 0.3); after
40 members, the performance of the EnKF hardly improvesas the number of members increases, and in the experimentswith 300 members, pure EnKF shows the worst score. Here,
since the sampling error decreases as the number of particles
increases, it would seem that using more of the local PF update would give better results. However, for this particularmodel and a linear observation operator, this is not the case.
Even with Ne 5 300, the best performing experiments use afactorization that amounts to 70% of the EnKF increment
being used. The experiment that determines the value ofk adaptively shows less optimal but suitable results for smaller
ensemble sizes}with the added beneﬁt of not needing to be
tuned. As the number of members increases, however, the sample size for SWT increases, thus making the test more accurate.
Increasing the ensemble size also increases the rejection rateof the null hypothesis, which is a desirable property. While themean value of k becomes larger as the ensemble size increases,
the value converges slowly to 0.2 (Fig. 8a). The L96 priors remainclose to Gaussian for most data assimilation cycles when using a
sufﬁciently dense network of observations with linear measurement operators, thus leading us to conclude that SWT operates
appropriately for this application.

F3

JOBNAME: MWR 00#00 2022 PAGE: 11 SESS: 8 OUTPUT: Sat Dec 24 08:35:30 2022 Total No. of Pages: 21/ams/mwr/0/mwrD220108

Monthly Weather Review (Proof Only)

MONTH 2022

KUROSAWA AND POTERJOY

11

FIG. 4. As in Fig. 3, but for (a) Mix Case1, (b) Mix Case2, and (c) Mix Case3.

Results from Nonlinear Case 1 are shown in Fig. 3b. With
Ne 5 20–40, the trend is the same as in the Linear Case: theexperiments with a more signiﬁcant percentage of the EnKF
show better scores. However, this feature is maintained evenwith large ensemble sizes. This is because of the precision andfrequency of the observations compared to ones in Linear
Case, as described in Kurosawa and Poterjoy (2021). Sincemodel variables are around the magnitude of O(10), the non^ +x^ with sy 5 1.0 provides very prelinear operator H(x) 5 x
cise information to characterize the posterior estimation. Thisfact, combined with the frequency of measurements, makes
Gaussian estimation more appropriate, as forecasts yield priormembers that are generally close to the truth. Therefore, we
can conﬁrm that k in the adaptive experiment uses a larger
percentage of the EnKF than in Linear Case (see Fig. 8a).
Figure 3c shows the mean RMSEs from experiments thatuse measurements simulated with Nonlinear Case 2. In experiments using this observation network, a situation occurs
in which the nonlinearity in the application becomes muchlarger than the sampling error in the prior and posterior distributions estimated by the ensemble. Owing to the strong nonlinearity of the observations, experiments using mainly theGaussian-based method struggle to provide an accurate
RMSE. In particular, the pure EnKF diverges, even with
Ne 5 300. The mean value of k in the adaptive experimentshows that most of the update is used for the local PF (Fig. 8a).
This result occurs as the strongly nonlinear observation operator tends to induce skewness in prior distributions, and SWT
frequently rejects the null hypothesis.
Based on the above results from the “normal” setting, resultsfrom “mixed” observation networks yield intuitive results

Fig(s). 4 live 4/C(Fig. 4). For example, in the case where the observation oper^ and H(x) 5 x^+x^ are combined, the experiators H(x) 5 x
ments with larger values of k tend to produce worse scores(Fig. 4a). On the other hand, experiments using strongly nonlinear operators experienced the best performance with smallvalues of k, and became unstable when the EnKF contribution was too large. As such, we consistently ﬁnd that the experiments with a value of k close to 0.5 are very stable. The
partial update by the local PF adjusts particles to a Gaussianlike distribution, providing an optimal prior distribution for
the EnKF update. The adaptive experiment also shows asatisfactory performance in the case of any combination ofobservation operators. We can see that using SWT is able
to estimate the optimal k according to each observation
operator.
We summarize the results for the “data void” setting in
Fig. 5. Mean RMSEs are uniformly higher than the “normal”
setting with any observation operators despite using thesmaller localization scale. Notably, in Nonlinear Case 2, several ﬁxed experiments diverged, but the experiments with theappropriate blend of the EnKF and the local PF (k 5 0.3–1.0)
are stable (Fig. 5c). The mean value of k used in the adaptive
experiment is close to 1 (Fig. 8b). We note that this experiment shows slight advantages over experiments that keep kﬁxed near 1, thus underscoring the importance of allowing k
to change over space and time.
Last, results obtained from simulated “model error” experiments show elevated errors for all experiments, regardless ofobservation operators in Fig. 6. The presence of model errors
means that the prior variance can be quite large, which leads
to more frequent non-Gaussian prior distributions for L96.

F4

F5

F6

JOBNAME: MWR 00#00 2022 PAGE: 12 SESS: 8 OUTPUT: Sat Dec 24 08:35:30 2022 Total No. of Pages: 21/ams/mwr/0/mwrD220108

Monthly Weather Review (Proof Only)

12

MONTHLY WEATHER REVIEWVOLUME 00

FIG. 5. As in Fig. 3, but using the data void observation network.

The hybrid strategies with a value of k close to 0.5 show
clear advantages in this regime. The parametric (Gaussian)assumption that follows the PF steps in hybrid conﬁgurations allows the ﬁlter to more easily adjust solutions for observations that lie outside the span of the ensemble. Hence,

it shifts particles closer to observations in a manner that is
not permitted by the PF}for variables that are detected tohave Gaussian errors.
To investigate the behavior of SWT speciﬁcation of k forthese simulations, we examine a sample time series of prior

FIG. 6. As in Fig. 3, but using an imperfect L96 model for forecast steps.

Fig(s). 5,6 live 4/C

JOBNAME: MWR 00#00 2022 PAGE: 13 SESS: 8 OUTPUT: Sat Dec 24 08:35:30 2022 Total No. of Pages: 21/ams/mwr/0/mwrD220108

MONTH 2022

KUROSAWA AND POTERJOY

13

FIG. 7. Time series of ensemble spread (red) and estimated k (blue) in adaptive hybrid experiments for Linear Case with “model error” setting. Values are from the ﬁrst variable of the L96
model in the ﬁrst trial of the experiments with ensemble size Ne 5 300. The correlation coefﬁcient between the pair of time series is indicated in top right.

Monthly Weather Review (Proof Only)

F7

F8

ensemble variance and estimated k for the experiment using alinear measurement operator (Fig. 7). The plotted values
come from the ﬁrst variable of the L96 model in the ﬁrst trialof the experiments with ensemble size Ne 5 300. Because of
an imperfect model, the prior variance ﬂuctuates signiﬁcantly
over the entire period. When the prior distribution has a larger
variance, the nonlinear model dynamics can more readily produce non-Gaussian priors, which SWT successfully detects. For
this example, the only factor that can contribute to non-Gaussianpriors is the nonlinear model itself as the measurement operator
is linear. Hence, ﬂuctuations of the ensemble variance and k arehighly correlated. In Fig. 7, the time- and space-average value ofk over the period is 0.2655, which is very close to the value
0.2766 in Fig. 8b (note that the value of k in Fig. 7 is from the ﬁrst
variable in the L96 model, while in Fig. 8 is from the average ofall variables in the model). As in the “data void” simulations, theexperiments with adaptively estimated k again show improvements over experiments with values of k that are conﬁgured touse close to the average mean estimated k, but ﬁxed over space
and time (Fig. 6a). In general, we ﬁnd that choosing k adaptively
is beneﬁcial in “model error” experiments, owing to its ability tomaintain ﬁlter stability without rigorous tuning. The sporadic
non-Gaussian priors produced by L96 in “model error” experiments introduce a major challenge that mimics the expected behavior of real weather systems.
Based on the above results, the statistical hypothesis testing
approach yields adequate hybrid factor estimates in all situations we examined for this study. Moreover, the approach hassigniﬁcant value for more realistic applications, such as nonhomogeneous observation networks and unknown model process error. We expect similar beneﬁts for geophysical problems
that are characterized by a variety of dynamic instabilities aswell. Furthermore, the proposed adaptive hybrid method avoids
the need to tune heuristic parameters, such as the hybrid factor,
which we ﬁnd to be sensitive to observation operators, observation density, and model process uncertainty.

FIG. 8. Mean estimated k of the adaptive hybrid experiments as a function of ensemble size.
Results are shown for (a) “normal” and “mix” settings in Figs. 3 and 4 and (b) “data void”
and “model error” settings in Figs. 5 and 6. Values are from the average of the last 9 years with
100 parallel trials.

Fig(s). 7,8 live 4/C

JOBNAME: MWR 00#00 2022 PAGE: 14 SESS: 8 OUTPUT: Sat Dec 24 08:35:31 2022 Total No. of Pages: 21/ams/mwr/0/mwrD220108

14

MONTHLY WEATHER REVIEWMonthly Weather Review (Proof Only)c. Idealized vortex model
In contrast to the low-dimensional applications used in theprevious subsections, realistic atmospheric forecast models
have several variables at each grid point, such as air temperature, winds, pressure, and speciﬁc humidity. These variables
also exhibit large spatial error dependence with one another,which is not accounted for in adaptive choices for k. As such,an observation of one variable must be used to update all collocated and nearby variables; for EnKFs, this step considers
prior error covariance across each variable. Furthermore, extending k to be the same dimension as the full state vector,
rather than the grid dimension, would bring additional algorithmic complexity to the proposed hybrid ﬁlter. To address
the problem of collocated variables for estimating k, a natural
choice is to perform the hypothesis test using all variables thatare expected to be correlated with that variable at a grid
point, i.e., by performing a test for multivariate normality.
The numerical experiments performed in this section serve
the purpose of illustrating the advantages of optimally adjusted k estimated via SWT extended to the test for multivariate normality proposed by Royston (1983). Considering
marginal PDFs in the test for multivariate normality is expected to provide a reasonable k for data assimilation updates
that account for correlations between collocated variables,
which is an important practical feature of the proposed hybrid
method in real weather applications, as the transition between
the local PF and EnKF updates are decided across grid points
rather than state variables. The method, however, still neglects dependence for variables at different grid points, which
is one theoretical shortcoming.
The current section provides an illustrative comparison of
the EnKF, the local PF, and hybrid updates using a low dimensional application that mimics a common challenge for ﬁltering geophysical ﬂow, namely, the problem of assimilating
measurements for mesoscale weather features that are notwell constrained by measurements.

1) EXPERIMENTAL DESIGNS
Adopting the same application introduced in Poterjoy (2022a),
we will reproduce the data assimilation challenge posed byalignment errors associated with mesoscale weather featuresby modeling a vortex wind ﬁeld with a Rankine vortex proﬁle(Acheson 1990). The model produces a case of a single vortex inzero-mean ﬂow, but position uncertainty. The Rankine vortex
consists of a wind ﬁeld exhibiting uniform vorticity in the vortex,
and an outer region of zero vorticity. For cylindrical coordinateswith the origin chosen to be the vortex center, and assuming that
all nonzero vorticity is uniformly distributed within a circle ofradius Rmax, the tangential winds uu are a function of radius r:U r , u Rmax
uu 5 
R Uu max ,rr , Rmax,(14)

r $ Rmax

where Uu is maximum wind speed. Both the radial wind component (ur) and vertical wind component (uz) are assumed to be

VOLUME 00

zero. For this demonstration, we transform winds into Cartesian
coordinates so that the model state vector comprises zonal (u)and meridional (y) wind components; i.e., x 5 (u, v)T.
The current study generates vortices on a two-dimensional
Cartesian region consisting of 91 3 91 equally spaced grid
points. To generate a prior sample, we ﬁrst designate a control
5 46,
state with the center of the vortex located at (iCTRLc5 46), and UuCTRL 5 30 m s21 and RCTRLjCTRLcmax 5 12. Then,the position and wind parameters of each vortex are randomly
drawn independently from a Gaussian distribution, and added
to the control vortex parameters. That is, the center of each, s2p ) and N(jCTRL, s2p ),
vortex (inc , jnc ) is sampled from N(iCTRLccrespectively, for n 5 1, … , Ne, where sp is a prescribed
position error standard deviation that changes for each prior.
Uun and Rnmax of each vortex is drawn from N(UuCTRL , 1) andN(RCTRL
max , 1), respectively.
Observations are generated by uniformly selecting points
from within the scan area of a hypothetical Doppler radarwith a radius of 30 grids, placed at coordinates (iradar 5 25,
jradar 5 25) in the lower-left corner of the domain. We produce each observation by projecting the truth state wind in
the direction of the hypothetical radar beam pointing outward
from the radar. In this experiment, the errors added to eachobservation are drawn from N(0, s2o ) for so 5 3, and we set
the number of observations Ny to be 100. Figure 9a shows thevalue of uu for the cross section through the center of the control state. Reproduced from Poterjoy (2022a), Fig. 9b shows asingle 15 m s21 wind speed contour for the Rankine vortex on
the 2D domain. The scan region by the virtual radar is indicated by the curved segment in the lower left part of the domain, which covers only one quadrant of the vortex, and thegreen and red dots indicate the location and magnitude of themeasurements (Fig. 9b).
As in the experiment with the L96 in section 4b, out of anabundance of caution, we perform 3000 parallel trials withunique sets of priors, true solutions, and observations in orderto capture the range of plausible outcomes for this application. To perform each trial, for the truth state, we generate
the center (itc , jtc ), maximum wind speed Uut , and radius Rt
from a Gaussian distribution as well as other prior ensemble
members. This process allows us to create a state in which thetrue value is indistinguishable from any prior member, i.e.,
the true value is also a sample from the prior distribution withequal probability, which is a condition assumed when performing data assimilation for real atmospheric models. The
truth state, which varies by each trial, is used to generate the observations and provides a reference for evaluating data assimilation experiments performed for each trial. We repeat thesetrials for several choices of position error standard deviation
sp 5 {0.0, 4.0, 8.0, 12.0} and ensemble size Ne 5 {40, 100, 300}.
For reference, Fig. 10 shows the variability of the initial ensemble members according to each sp with Ne 5 40.
We performed a total of 12 experiments, one in which thevalue of k is estimated adaptively, and the others in which k is
spatially constant between 0 and 1, using increments of 0.1. In
the adaptive hybrid experiment, we use SWT extended to thetest for multivariate normality; that is, the Gaussianity in theprior samples is detected using two variables, u and y. In this

F9

F10

JOBNAME: MWR 00#00 2022 PAGE: 15 SESS: 8 OUTPUT: Sat Dec 24 08:35:31 2022 Total No. of Pages: 21/ams/mwr/0/mwrD220108

KUROSAWA AND POTERJOYMonthly Weather Review (Proof Only)

MONTH 2022

FIG. 9. (a) Tangential wind speed as a function of grid points, calculated using the Rankine vortex model with point 46 as the center
location. (b) The 15 m s21 wind speed contours for the vortex
placed on a 2D grid; values greater than 15 m s21 are indicated by
the hatched region. The green and red markers indicate the location and magnitude of radial wind observations created for a synthetic radar located at coordinate (iradar, jradar). This ﬁgure is a reproduction of Fig. 8 in Poterjoy (2022a).

demonstration, we assimilate the observations in each of the ksettings using each of the choices of prior and then calculate
the RMSEs of the posterior mean relative to the true wind velocity. All experiments use the localization function f:
1
22 



 1 d(i,j) ,
f 5 exp2



 2 s


(15)

loc

where d(i, j) is the physical distance between grid points i and j,and sloc is the localization parameter to scale the width of localization, which is set to 2000 in the current study.

2) RESULTS
This subsection discusses results obtained by performing
single-cycle data assimilation using the described sets ofFig(s). 9 live 4/C

15

observations and prior members. Figure 11 shows the posterior
RMSEs for four experiments as examples: EnKF (k 5 0.0),
the local PF (k 5 1.0), PF–EnKF with k 5 0.5, and PF–EnKF
with adaptive k estimation. For experiments that use sp 5 0.0
for the prior, each data assimilation method shows low
RMSEs that are visibly similar (ﬁrst row of Fig. 11). Thisﬁnding is expected because the small sp 5 0.0 leads to theGaussian assumption being valid (Poterjoy 2022a). However,as the value of sp increases, each experiment yields vast differences in the upper right portion of the domain where eachﬁlter must infer wind estimates from distant observations.
This application is especially problematic for the EnKF, as linear updates do not properly capture nonlinear dependence inwinds across the vortex (Poterjoy 2022a). While the local PF
and PF–EnKF produce smaller mean RMSEs than the EnKF,
we note that these errors also continue to decrease as thenumber of members increases, because of the decrease insampling error. Furthermore, compared to the EnKF and local PF experiments, both hybrid experiments show more accurate results across all domains, which demonstrates that the
hybrid PF–EnKF method is effective at shifting particles into
an approximate Gaussian before applying the EnKF step}even for the highly non-Gaussian vortex application discussed
in Poterjoy (2022a).
Comparing experiments with ﬁxed and adaptive k, Fig. 12shows the grid points over the domain where a speciﬁed valueof k produced the smallest RMSEs for each choice of prior.
For all settings, the adaptive estimate yields the smallest
errors outside of the vortex, thus reﬂecting diversity in optimal k in this region. When sp 5 0, experiments with a valueof k ﬁxed near 0.5 show the best results near the vortex center(Figs. 12a–c). We suspect this result occurs because the location of the prior and true vortex centers are the same for all
prior distributions when sp 5 0, but since R is drawn from
N(RCTRL, 1), the winds exhibit bimodal behavior, which iscontrolled by parameter R in (14); recall, this parameter
divides the domain into regions of zero and nonzero}but
constant}vorticity. Since the region near the vortex center ischaracterized by the presence of both zeros and nonzeros, it
is conceivable that the case of k 5 0.5, where both PF and
EnKF can be used in a balanced manner, happens to be themost optimal. Therefore, a ﬁxed value of k can be identiﬁed
via rigorous tuning, rather than resorting to hypothesis testing. When sp . 0 and a sufﬁciently large ensemble size isused (e.g., Ne 5 300), SWT correctly identiﬁes values for kthat outperform ﬁxed values for k over most of the domain(Figs. 12f,i,l). Prior vortices are no longer collocated as sp increases, so the region where ﬁxed values for k are optimal
gradually extends outward from the center.
We also examine the mean value of k (averaged over trials)
in the experiments where k is adaptively adjusted (Fig. 13).
First, for the case of sp 5 0.0, the area close to the center ofthe control vortex, where experiments with a value of k ﬁxed
near 0.5 corresponds to the location where estimated k isabout 0.2–0.5 in the ﬁrst row of Fig. 13. This indicates that thearea has a low percentage of the local PF updates comparedto the ﬁxed experiment. For the cases of Ne 5 40 and 100, theareas where the experiment using the adaptively adjusted k is

F11

F12

F13

JOBNAME: MWR 00#00 2022 PAGE: 16 SESS: 8 OUTPUT: Sat Dec 24 08:35:31 2022 Total No. of Pages: 21/ams/mwr/0/mwrD220108

Monthly Weather Review (Proof Only)

16

MONTHLY WEATHER REVIEWVOLUME 00

FIG. 10. Variability of the initial ensemble members for (a) sp 5 0.0, (b) sp 5 4.0, (c) sp 5 8.0, and (d) sp 5 12.0 with ensemble size
Ne 5 40. Each colored line shows the 15 m s21 wind speed contours. The black dashed lines show the 15 m s21 wind speed contours of thecontrol state.

inferior in Fig. 12 are generally estimated to have a value of kless than 0.5 (the ﬁrst two columns of Fig. 13). However, inthe case with Ne 5 300, the values of k in those locations are
generally more than 0.5, and the difference from the ﬁxed experiment is not signiﬁcant (third column of Fig. 13). This may
be because the larger sample size used in SWT leads to more
frequent rejection of the null hypothesis, resulting in moreiterations of the local PF. Furthermore, the experiment withestimated k is generally more stable in the other areas, faraway from the center, especially in the upper-right domains,
where there are no observations.
As in the previous subsections, the kinematic vortex experiments illustrate the advantage of specifying k adaptively versus keeping this parameter ﬁxed. In this demonstration,
however, k is estimated using u and y with SWT extended tohandle multivariate normality. Hence, results from this section differ in that we successfully analyzed samples from multivariate probability distributions. The multivariate approach
is more practical for data assimilation with real atmosphericmodels, when we need to consider correlation among collocated variables. Furthermore, as noted by Poterjoy (2022a),the non-Gaussian data assimilation problem constructed inthis section has univariate marginal distributions that areclose to Gaussian, but multivariate marginals for variables
across grid points that are far from Gaussian. In terms ofcomputing burden, high-resolution models, such as thoseused for weather forecasting, are constrained by ensemble
size. As a result, detecting such characteristics with a limited number of ensembles to be utilized for operationalmodels is extremely challenging. Nevertheless, we ﬁnd the
proposed SWT approach to be sufﬁcient for identifying deviations from a multivariate normality for collocated winds,
which shows added value over a rigorously tuned hybrid
methodology that uses ﬁxed speciﬁcations for k.

5. Discussion and conclusions
The current study introduces a novel approach to forming
an adaptive hybrid data assimilation method that mixes thetheoretical strengths and ﬂexibility of particle ﬁlters with

Fig(s). 10 live 4/C

Gaussian-based ensemble Kalman ﬁlters (EnKFs), which are
more resilient to bias in sample-estimated prior uncertainty.
For this purpose, we use a recently proposed PF by Poterjoy(2022b), which introduces a regularization and tempering
methodology to improve ﬁlter performance when sampling
error is large. The tempering step consists of a factorization of
the particle weights, which provides a natural framework for
combining local PFs with alternative ﬁlters to mitigate the effects of sampling error. In addition to identifying portions of
the state space where a PF may provide more accurate marginal posterior estimates than an EnKF, the adaptive strategy
can switch between ﬁlters partway through data assimilation
steps. The latter property is beneﬁcial when Gaussian assumptions are appropriate for posteriors but not for priors, whichis common when likelihoods are Gaussian. In this case, partial
updates performed by the PF can adjust the distribution ofparticles to more closely ﬁt a Gaussian, which allows for amore effective use of EnKFs. To determine the timingof the transition between these ﬁlter updates, we use theShapiro–Wilk test (SWT), which has excellent power among
omnibus tests to detect deviations from normality. The useof SWT allows for accurate detection of Gaussianity evenwhen the ensemble size is small. SWT also requires minimal
computing time, thus permitting its use between PF iterations, which can be carried out until prior sample distributions for marginals at each grid point are detected to be
Gaussian. Increasing the ensemble size also increases therejection rate of the null hypothesis and leads to a smaller
portion of updates being made by an EnKF, which is a
desirable property.
To examine the performance of the adaptive hybrid, thisstudy constructs numerous data assimilation experiments using
a low-dimensional dynamical model, which is characterized by
40 equally spaced variables on a periodic domain. In general,
the statistical hypothesis testing approach yields adequate estimates of the hybrid factor in all situations considered in thisstudy. Given a homogeneous network of equally spaced observations, the adaptive formulations of the hybrid ﬁlter are as accurate as the rigorously tuned hybrid parameters. The adaptive

JOBNAME: MWR 00#00 2022 PAGE: 17 SESS: 8 OUTPUT: Sat Dec 24 08:35:31 2022 Total No. of Pages: 21/ams/mwr/0/mwrD220108

KUROSAWA AND POTERJOY

17Monthly Weather Review (Proof Only)

MONTH 2022

FIG. 11. Analysis RMSEs of velocity for 1) EnKF (k 5 0.0), 2) PF (k 5 1.0), 3) the experiment with k 5 0.5, and 4) with adaptive estimated k.
The position error standard deviation sp is (a)–(c) 0.0, (d)–(f) 4.0, (g)–(i) 8.0, and (j)–(l) 12.0. The ensemble size Ne is (left) 40; (center) 100; and(right) 300. Values are from the average of 3000 parallel trials. The black dashed lines show the 15 m s21 wind speed contours of the control state.

approach also demonstrates clear advantages in experiments
containing heterogeneous observation networks and unknown
model process errors}in which case, the optimal choice of adaptive parameter varies temporally or across variables.

Fig(s). 11 live 4/C

The study also examines practical challenges for adopting
the new method for real Earth system models, which are characterized by multiple variables at common grid points andlarge error correlations through space; e.g., modern weather

JOBNAME: MWR 00#00 2022 PAGE: 18 SESS: 8 OUTPUT: Sat Dec 24 08:35:32 2022 Total No. of Pages: 21/ams/mwr/0/mwrD220108

MONTHLY WEATHER REVIEWVOLUME 00

Monthly Weather Review (Proof Only)

18

FIG. 12. The experiment with the lowest RMSEs out of the 12 different k cases over 3000 parallel trials. Color represents the best kcase, which produced the smallest RMSE, comparing ﬁxed and adaptive k experiments. The black dashed lines show the 15 m s21 windspeed contours of the control state.

prediction models. The computational expense of such models limits the amount of tuning that can be performed for heuristic parameters used during data assimilation, which can besensitive to observation operators, observation frequency, andmodel process uncertainty. Therefore, this study adopts an idealized kinematic vortex model to study the behavior of theadaptive hybrid. This model permits large error dependenceFig(s). 12 live 4/C

across variables displaced over a two-dimensional domain,and contains two variables (zonal and meridional winds) ateach grid point thus requiring a multivariate SWT to adaptively choose how to partition PF and EnKF updates. For thisapplication, the hybrid factor is estimated using SWT extended
to detect multivariate normality for ensembles of u and yat each grid point. This approach allows the use of the

JOBNAME: MWR 00#00 2022 PAGE: 19 SESS: 8 OUTPUT: Sat Dec 24 08:35:32 2022 Total No. of Pages: 21/ams/mwr/0/mwrD220108

KUROSAWA AND POTERJOY

19Monthly Weather Review (Proof Only)

MONTH 2022

FIG. 13. Mean value of k in the experiments with adaptive estimated k. Values are from the average of 3000 parallel trials. The black
dashed lines show the 15 m s21 wind speed contours of the control state.

appropriate factor to account for multivariate marginal
distributions for updating the state variables, alongsideobservation-space priors. Specifying the hybrid factor for collocated variables also simpliﬁes the algorithmic formulationof the adaptive methodology, as it only requires the factor to
be speciﬁed for all grid points and observation-space priors
used during data assimilation. The experiments reveal spatial

Fig(s). 13 live 4/C

patterns of adaptively chosen hybrid factors that result inlarge PF updates in portions of the state space where Gaussian assumptions are known to be incorrect, and are close
to the values identiﬁed at each grid point from rigorously
tuned experiments aimed at reducing posterior mean RMSEs.
These results encourage further testing for real geophysical problems that are characterized by a variety of dynamic instabilities.

JOBNAME: MWR 00#00 2022 PAGE: 20 SESS: 8 OUTPUT: Sat Dec 24 08:35:33 2022 Total No. of Pages: 21/ams/mwr/0/mwrD220108

Monthly Weather Review (Proof Only)

20

MONTHLY WEATHER REVIEWIn summary, the proposed adaptive hybrid method performs
well in idealized simulations that mimic data assimilation problems encountered for real geophysical modeling systems. Because
the new strategy relies on statistical hypothesis testing, it becomes
more stable when the ensemble size increases. The proposed
method obviates the need for tuning a hybrid parameter that inﬂuences when an EnKF is preferred over PF, which can depend
on a number of factors including the underlying model dynamics
and observation network. This property of the method has theoretical beneﬁts for real Earth system models where rigorous tuning of data assimilation parameters is not always feasible, and theshape of error distributions is ﬂow dependent. Last, this study
demonstrates how SWT can be extended to consider error dependence for collocated variables. Further research will explore
the use of multivariate error dependence for variables across grid
points, which may be needed for prior distributions that are characterized by strong nonlinear dependence for variables displaced
geographically.
Acknowledgments. Funding for this work was provided
by NOAA Grant NA20OAR4600281 and NSF/CAREER
Award AGS1848363.
Data availability statement. This paper uses numerous lowdimensional model simulations that can be easily replicatedusing software written for this study. The model code, compilation script, and the namelist settings are available at https://github.com/Kenta9638/MWR_2022.
REFERENCES
Acheson, D. J., 1990: Elementary Fluid Dynamics. Oxford University Press, 408 pp.
Althouse, L. A., W. B. Ware, and J. M. Ferron, 1998: Detecting
departures from normality: A Monte Carlo simulation of a newomnibus test based on moments. Diversity and Citizenship in
Multicultural Societies, San Diego, CA, American Educational
Research Association, 33 pp., https://eric.ed.gov/?id=ED422385.
Anderson, T. W., and D. A. Darling, 1954: A test of goodness ofﬁt. J. Amer. Stat. Assoc., 49, 765–769, https://doi.org/10.1080/
01621459.1954.10501232.
Arshad, M., M. T. Rasool, and M. I. Ahmad, 2003: Anderson
Darling and modiﬁed Anderson Darling tests for generalized
Pareto distribution. J. Appl. Sci., 3, 85–88, https://doi.org/10.
3923/jas.2003.85.88.
Bocquet, M., C. A. Pires, and L. Wu, 2010: Beyond Gaussianstatistical modeling in geophysical data assimilation. Mon. Wea.
Rev., 138, 2997–3023, https://doi.org/10.1175/2010MWR3164.1.
Chustagulprom, N., S. Reich, and M. Reinhardt, 2016: A hybrid
ensemble transform particle ﬁlter for nonlinear and spatiallyextended dynamical systems. SIAM/ASA J. Uncertainty
Quantif., 4, 592–608, https://doi.org/10.1137/15M1040967.
Cramér, H., 1928: On the composition of elementary errors.
Scand. Actuar. J., 1928, 13–74, https://doi.org/10.1080/03461238.
1928.10416862.
Dufour, J., A. Farhat, L. Gardiol, and L. Khalaf, 1998: Simulationbased ﬁnite sample normality tests in linear regressions. Econom.
J., 1, C154–C173, https://doi.org/10.1111/1368-423X.11009.
Emerick, A. A., and A. C. Reynolds, 2012: History matchingtime-lapse seismic data using the ensemble Kalman ﬁlter withVOLUME 00

multiple data assimilations. Comput. Geosci., 16, 639–659,
https://doi.org/10.1007/s10596-012-9275-5.
Errico, R. M., P. Bauer, and J.-F. Mahfouf, 2007: Issues regarding
the assimilation of cloud and precipitation data. J. Atmos.
Sci., 64, 3785–3798, https://doi.org/10.1175/2006JAS2044.1.
Evensen, G., 1994: Sequential data assimilation with a nonlinear
quasi-geostrophic model using Monte Carlo methods toforecast error statistics. J. Geophys. Res., 99, 10 143–10 162,https://doi.org/10.1029/94JC00572.
}}, and P. J. van Leeuwen, 2000: An ensemble Kalmansmoother for non-linear dynamics. Mon. Wea. Rev., 128,
1852–1867, https://doi.org/10.1175/1520-0493(2000)128,1852:
AEKSFN.2.0.CO;2.
Fabry, F., and J. Sun, 2010: For how long should what data be assimilated for the mesoscale forecasting of convection andwhy? Part I: On the propagation of initial condition errors
and their implications for data assimilation. Mon. Wea. Rev.,
138, 242–255, https://doi.org/10.1175/2009MWR2883.1.
Farrell, P. J., M. Salibian-Barrera, and K. Naczk, 2007: On tests
for multivariate normality and associated simulation studies.
J. Stat. Comput. Simul., 77, 1065–1080, https://doi.org/10.1080/
10629360600878449.
Frei, M., and H. R. Kunsch, 2013: Bridging the ensemble Kalman
and particle ﬁlters. Biometrika, 100, 781–800, https://doi.org/
10.1093/biomet/ast020.
Gaspari, G., and S. E. Cohn, 1999: Construction of correlation
functions in two and three dimensions. Quart. J. Roy. Meteor.
Soc., 125, 723–757, https://doi.org/10.1002/qj.49712555417.
Grooms, I. G., and G. Robinson, 2021: A hybrid particle-ensemble
Kalman ﬁlter for problems with medium nonlinearity. PLOS
ONE, 16, e0248266, https://doi.org/10.1371/journal.pone.0248266.
Honda, T., and Coauthors, 2018: Assimilating all-sky Himawari-8
satellite infrared radiances: A case of Typhoon Soudelor(2015). Mon. Wea. Rev., 146, 213–229, https://doi.org/10.1175/
MWR-D-16-0357.1.
Houtekamer, P. L., and H. L. Mitchell, 1998: Data assimilation using
an ensemble Kalman ﬁlter technique. Mon. Wea. Rev., 126,
796–811, https://doi.org/10.1175/1520-0493(1998)126,0796:
DAUAEK.2.0.CO;2.
Kolmogorov, A. N., 1933: Sulla determinazione empirica di una
legge di distribuzione. Giorn. Instit. Ital. Attu., 4, 83–91.
Kondo, K., and T. Miyoshi, 2019: Non-Gaussian statistics in global
atmospheric dynamics: A study with a 10 240-member ensemble Kalman ﬁlter using an intermediate atmospheric general circulation model. Nonlinear Processes Geophys., 26,
211–225, https://doi.org/10.5194/npg-26-211-2019.
Kullback, S., and R. A. Leibler, 1951: On information and sufﬁciency. Ann. Math. Stat., 22, 79–86, https://doi.org/10.1214/
aoms/1177729694.
Kurosawa, K., and J. Poterjoy, 2021: Data assimilation challenges
posed by nonlinear operators: A comparative study of ensemble and variational ﬁlters and smoothers. Mon. Wea. Rev., 149,
2369–2389, https://doi.org/10.1175/MWR-D-20-0368.1.
Li, R., N. Magbool Jan, B. Huang, and V. Prasad, 2019:
Constrained multimodal ensemble Kalman ﬁlter based on
Kullback–Leibler (Kl) divergence. J. Process Control, 79,
16–28, https://doi.org/10.1016/j.jprocont.2019.03.012.
Lilliefors, H. W., 1967: On the Kolmogorov-Smirnov test for normality with mean and variance unknown. J. Amer. Stat. Assoc.,
62, 399–402, https://doi.org/10.1080/01621459.1967.10482916.
Lorenz, E. N., 1996: Predictability: A problem partly solved. Proc.
Seminar on Predictability, Shinﬁeld Park, Reading, United
Kingdom, ECMWF, 1–18, https://www.ecmwf.int/node/10829.

JOBNAME: MWR 00#00 2022 PAGE: 21 SESS: 8 OUTPUT: Sat Dec 24 08:35:34 2022 Total No. of Pages: 21/ams/mwr/0/mwrD220108

Monthly Weather Review (Proof Only)

MONTH 2022

KUROSAWA AND POTERJOY}}, and K. A. Emanuel, 1998: Optimal sites for supplementary
weather observations: Simulation with a small model. J. Atmos.
Sci., 55, 399–414, https://doi.org/10.1175/1520-0469(1998)055,0399:OSFSWO.2.0.CO;2.
Mendes, M., and A. Pala, 2003: Type I error rate and power ofthree normality tests. Inf. Technol. J., 2, 135–139, https://doi.
org/10.3923/itj.2003.135.139.
Metref, S., E. Cosme, C. Snyder, and P. Brasseur, 2014: A nonGaussian analysis scheme using rank histograms for ensemble
data assimilation. Nonlinear Processes Geophys., 21, 869–885,
https://doi.org/10.5194/npg-21-869-2014.
Minamide, M., and F. Zhang, 2017: Adaptive observation error
inﬂation for assimilating all-sky satellite radiance. Mon. Wea.
Rev., 145, 1063–1081, https://doi.org/10.1175/MWR-D-16-0257.1.
Morzfeld, M., and D. Hodyss, 2019: Gaussian approximations inﬁlters and smoothers for data assimilation. Tellus, 71A, 1–27,
https://doi.org/10.1080/16000870.2019.1600344.
}}, }}, and J. Poterjoy, 2018: Variational particle smoothers
and their localization. Quart. J. Roy. Meteor. Soc., 144, 806–
825, https://doi.org/10.1002/qj.3256.
Nerger, L., 2022: Data assimilation for nonlinear systems with ahybrid nonlinear Kalman ensemble transform ﬁlter. Quart. J.
Roy. Meteor. Soc., 148, 620–640, https://doi.org/10.1002/qj.4221.
Penny, S. G., and T. Miyoshi, 2016: A local particle ﬁlter for highdimensional geophysical systems. Nonlinear Processes Geophys., 23, 391–405, https://doi.org/10.5194/npg-23-391-2016.
Pimentel, S., and Y. Qranfal, 2021: A data assimilation frameworkthat uses the Kullback-Leibler divergence. PLOS ONE, 16,e0256584, https://doi.org/10.1371/journal.pone.0256584.
Poterjoy, J., 2016: A localized particle ﬁlter for high-dimensional
nonlinear systems. Mon. Wea. Rev., 144, 59–76, https://doi.
org/10.1175/MWR-D-15-0163.1.
}}, 2022a: Implications of multivariate non-Gaussian data assimilation for multi-scale weather prediction. Mon. Wea. Rev., 150,
1475–1493, https://doi.org/10.1175/MWR-D-21-0228.1.
}}, 2022b: Regularization and tempering for a moment-matching
localized particle ﬁlter. Quart. J. Roy. Meteor. Soc., 148, 2631–
2651, https://doi.org/10.1002/qj.4328.
}}, and J. L. Anderson, 2016: Efﬁcient assimilation of simulated
observations in a high-dimensional geophysical system using
a localized particle ﬁlter. Mon. Wea. Rev., 144, 2007–2020,
https://doi.org/10.1175/MWR-D-15-0322.1.
}}, R. A. Sobash, and J. L. Anderson, 2017: Convective-scaledata assimilation for the Weather Research and Forecasting
Model using the local particle ﬁlter. Mon. Wea. Rev., 145,
1897–1918, https://doi.org/10.1175/MWR-D-16-0298.1.
}}, L. J. Wicker, and M. Buehner, 2019: Progress toward theapplication of a localized particle ﬁlter for numerical weather
prediction. Mon. Wea. Rev., 147, 1107–1126, https://doi.org/
10.1175/MWR-D-17-0344.1.
Potthast, R., A. Walter, and A. Rhodin, 2019: A localized adaptive
particle ﬁlter within an operational NWP framework. Mon. Wea.
Rev., 147, 345–362, https://doi.org/10.1175/MWR-D-18-0028.1.
Privé, N. C., Y. Xie, J. S. Woollen, S. E. Koch, R. Atlas, andR. E. Hood, 2013: Evaluation of the Earth Systems ResearchLaboratory’s Global Observing System Simulation Experiment
system. Tellus, 65A, 19011, https://doi.org/10.3402/tellusa.v65i0.
19011.
Razali, N., and Y. Wah, 2011: Power comparison of Shapiro–Wilk, Kolmogorov–Smirnov, Lilliefors and Anderson Darling
tests. J. Stat. Model. Anal., 2, 21–33.

21

Royston, J. P., 1982: An extension of Shapiro and Wilk’s W test
for normality to large samples. J. Roy. Stat. Soc., 31, 115–124,
https://doi.org/10.2307/2347973.
}}, 1983: Some techniques for assessing multivarate normality
based on the Shapiro–Wilk. J. Roy. Stat. Soc., 32, 121–133,
https://doi.org/10.2307/2347291.
}}, 1992: Approximating the Shapiro–Wilk W-test for nonnormality. Stat. Comput., 2, 117–119, https://doi.org/10.1007/
BF01891203.
}}, 1995: Remark AS R94: A remark on algorithm AS 181: The
W-test for normality. J. Roy. Stat. Soc., 44, 547–551, https://doi.
org/10.2307/2986146.
Ruiz, J., G.-Y. Lien, K. Kondo, S. Otsuka, and T. Miyoshi, 2021:
Reduced non-Gaussianity by 30 s rapid update in convectivescale numerical weather prediction. Nonlinear Processes Geophys., 28, 615–626, https://doi.org/10.5194/npg-28-615-2021.
Saculinggan, M., and E. A. Balase, 2013: Empirical power comparison of goodness of ﬁt tests for normality in the presenceof outliers. J. Phys. Conf. Ser., 435, 012041, https://doi.org/10.
1088/1742-6596/435/1/012041.
Shapiro, S. S., and M. B. Wilk, 1965: An analysis of variance test
for normality (complete samples). Biometrika, 52, 591–611,
https://doi.org/10.2307/2333709.
Slivinski, L., E. Spiller, A. Apte, and B. Sandstede, 2015: A hybrid
particle-ensemble Kalman ﬁlter for Lagrangian data assimilation. Mon. Wea. Rev., 143, 195–211, https://doi.org/10.1175/
MWR-D-14-00051.1.
Smirnov, N. V., 1939: Estimate of deviation between empirical
distribution functions in two independent samples. Bull. Moscow Univ., 2, 3616.
Srivastava, M., and T. Hui, 1987: On assessing multivariate normality based on Shapiro–Wilk W statistic. Stat. Probab. Lett.,
5, 15–18, https://doi.org/10.1016/0167-7152(87)90019-8.
Stengel, M., P. Unden, M. Lindskog, P. Dahlgren, N. Gustafsson,
and R. Bennartz, 2009: Assimilation of SEVIRI infrared radiances with HIRLAM 4D-Var. Quart. J. Roy. Meteor. Soc.,
135, 2100–2109, https://doi.org/10.1002/qj.501.
Stordal, A. S., H. A. Karlsen, G. Nævdal, H. J. Skaug, and B.
Vallès, 2011: Bridging the ensemble Kalman ﬁlter and particle
ﬁlters: The adaptive Gaussian mixture ﬁlter. Comput. Geosci.,
15, 293–305, https://doi.org/10.1007/s10596-010-9207-1.
Thadewald, T., and H. Büning, 2007: Jarque–Bera test and its competitors for testing normality}A power comparison. J. Appl.
Stat., 34, 87–105, https://doi.org/10.1080/02664760600994539.
Villaseñor Alva, J. A., and E. González-Estrada, 2009: A generalization of Shapiro–Wilk’s test for multivariate normality.
Commun. Stat. Theory Methods, 38, 1870–1883, https://doi.
org/10.1080/03610920802474465.
Vukicevic, T., T. Greenwald, M. Zupanski, D. Zupanski, T. V.
Haar, and A. S. Jones, 2004: Mesoscale cloud state estimationfrom visible and infrared satellite radiances. Mon. Wea. Rev.,
132, 3066–3077, https://doi.org/10.1175/MWR2837.1.
Whitaker, J. S., and T. M. Hamill, 2002: Ensemble data assimilation
without perturbed observations. Mon. Wea. Rev., 130, 1913–1924,
https://doi.org/10.1175/1520-0493(2002)130,1913:EDAWPO.2.0.
CO;2.
Zhang, F., C. Snyder, and J. Sun, 2004: Impacts of initial estimate
and observation availability on convective-scale data assimilation with an ensemble Kalman ﬁlter. Mon. Wea. Rev., 132,
1238–1253, https://doi.org/10.1175/1520-0493(2004)132,1238:
IOIEAO.2.0.CO;2.