JULY 2021

KUROSAWA AND POTERJOY

2369

Data Assimilation Challenges Posed by Nonlinear Operators: A Comparative Study of Ensemble
and Variational Filters and Smoothers
KENTA KUROSAWAa AND JONATHAN POTERJOYa,ba

Department of Atmospheric and Oceanic Science, University of Maryland, College Park, College Park, Marylandb
NOAA/Atlantic Oceanographic and Meteorological Laboratory, Miami, Florida(Manuscript received 6 November 2020, in final form 21 April 2021)

ABSTRACT: The ensemble Kalman filter (EnKF) and the 4D variational method (4DVar) are the most commonly usedfilters and smoothers in atmospheric science. These methods typically approximate prior densities using a Gaussian andsolve a linear system of equations for the posterior mean and covariance. Therefore, strongly nonlinear model dynamics andmeasurement operators can lead to bias in posterior estimates. To improve the performance in nonlinear regimes, minimization of the 4DVar cost function typically follows multiple sets of iterations, known as an ‘‘outer loop,’’ which helps
reduce bias caused by linear assumptions. Alternatively, ‘‘iterative ensemble methods’’ follow a similar strategy of periodically relinearizing model and measurement operators. These methods come with different, possibly more appropriate,assumptions for drawing samples from the posterior density, but have seen little attention in numerical weather prediction(NWP) communities. Last, particle filters (PFs) present a purely Bayesian filtering approach for state estimation, which
avoids many of the assumptions made by the above methods. Several strategies for applying localized PFs for NWP havebeen proposed very recently. The current study investigates intrinsic limitations of current data assimilation methodology
for applications that require nonlinear measurement operators. In doing so, it targets a specific problem that is relevant tothe assimilation of remotely sensed measurements, such as radar reflectivity and all-sky radiances, which pose challenges for
Gaussian-based data assimilation systems. This comparison includes multiple data assimilation approaches designed recently for nonlinear/non-Gaussian applications, as well as those currently used for NWP.
KEYWORDS: Kalman filters; Data assimilation; Nonlinear models; Numerical weather prediction/forecasting

1. Introduction
The ensemble Kalman filter (EnKF; Evensen 1994; Houtekamer
and Mitchell 1998; Evensen and van Leeuwen 2000) and thefour-dimensional variational method (4DVar; Thepáut andCourtier 1991) are the most commonly used filters andsmoothers in atmospheric science. Ensemble/variational hybrid approaches (e.g., Hamill and Snyder 2000; Lorenc 2003;
Buehner 2005) combine the flow-dependent ensemble covariance from an EnKF with climate-based covariance fromvariational methods. The methods have also become wellestablished and widely accepted for global weather prediction at major environmental prediction centers, such as theEuropean Centre for Medium-Range Weather Forecasts(ECMWF), the Met Office, the Environment and Climate
Change Canada (ECCC), and the National Centers forEnvironmental Prediction (NCEP). One strategy of the hybrid methods, denoted as ensemble-4DVar (E4DVar; Zhang
et al. 2009) in this manuscript, typically uses tangent linear
and adjoint model operators to minimize a cost function in
the same manner as the traditional 4DVar data assimilation
system. A second strategy is 4D-ensemble-Var (4DEnVar),in which the cost function minimization is computed based onan ensemble forecast instead of using tangent linear and adjoint models. In the 4DEnVar, temporal covariances are estimated from an ensemble of model trajectories that passthrough the observation time window. In either case, both

Corresponding author: Kenta Kurosawa, kkurosaw@umd.edu

methods approximate prior densities using a Gaussian andperform linearizations to relax these assumptions. Therefore,
strongly nonlinear model dynamics or measurement operators cause these methods to be biased, which leads to thesuboptimal use of major Earth observing systems, such as
satellite radiometers. For example, the combined impact ofhighly nonlinear model dynamics and measurement operators introduces major data assimilation challenges in weather
regimes containing clouds or precipitation. As a result, mostinfrared satellite assimilation studies mainly focus on clearsky observations (e.g., Errico et al. 2007; Fabry and Sun 2010;
Geer and Bauer 2011; Zou et al. 2013; Okamoto et al. 2014;
Minamide and Zhang 2017; Honda et al. 2018). This follows
despite the known benefits of assimilating cloudy radiances
for weather forecasting (e.g., Vukicevic et al. 2004; Stengel
et al. 2009; Privé et al. 2013). Some operational centers aremaking efforts to cope with these issues and assimilate cloudy
and precipitating microwave radiances (e.g., Zhu et al. 2016;
Geer et al. 2017, 2019). For further details on significant advances and current plans of operational centers that are close
to implementing assimilation, we encourage readers to review the summary presented in Geer et al. (2018).
Several procedures have been proposed to improve theperformance of these methods in nonlinear regimes. For example, in order to deal with issues within the 4DVar system(e.g., Bonavita et al. 2018), minimization of the 4DVar costfunction typically follows multiple sets of iterations to relinearize tangent linear and adjoints for the model, measurement operators, or both around an improved backgroundsolution. This step, known as an ‘‘outer loop,’’ helps reduce

DOI: 10.1175/MWR-D-20-0368.1
Ó 2021 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright
Policy (www.ametsoc.org/PUBSReuseLicenses).
Brought to you by University of Maryland, McKeldin Library | Unauthenticated | Downloaded 06/30/21 02:21 AM UTC

2370

MONTHLY WEATHER REVIEW

bias caused by linear assumptions, thus making Gaussian error
approximations more appropriate. The minimization strategy
follows the Gauss–Newton method, which is guaranteed toapproximate the posterior mode for local minima.
Alternatively, a number of methods fall under the genericcategory of ‘‘iterative ensemble methods,’’ which follow asimilar strategy of periodic relinearization. Note that here‘‘iterations’’ refers to multiple adjustments at a single time.
Both 4DVar and the iterative ensemble methods relinearize
the observation operator. The only difference is that in 4DVar,
the observation operator contains the nonlinear model. Gu andOliver (2007) introduced the ensemble randomized maximal
likelihood filter (EnRML) to handle nonlinearity by means of
iterations of the EnKF. Sakov et al. (2012) proposed the iterative ensemble Kalman filter (IEnKF), which uses a deterministic update form, ensemble square root filter, while
EnRML uses a stochastic update form, perturbed observations
method. Following the introduction of ensemble Kalmansmoother (EnKS; van Leeuwen and Evensen 1996; Evensen
and van Leeuwen 2000) for use in history matching bySkjervheim et al. (2011), the iterative forms of smoothers havedeveloped into useful tools by the reservoir-engineering community for history matching reservoir models. Chen and Oliver(2012) proposed an iterative form of EnRML targeted for oilreservoir modeling, and Bocquet and Sakov (2014) developed
the iterative ensemble Kalman smoother (IEnKS), which extends IEnKF using a fixed-lag smoother with an ensemble
variational method.
Emerick and Reynolds (2012) introduced the multiple data
assimilation scheme (MDA) to improve EnKF estimates fornonlinear cases by assimilating the same data multiple times
with the covariance matrix of the measurement errors multiplied by the number of data assimilation. We note that thename ‘‘MDA’’ is somewhat deceiving, as it is simply an application of tempering (Neal 1996). The process of the EnKF
with MDA (EnKF-MDA) is based on the idea that a ‘‘large
jump’’ between the forecast and analysis states could be reduced by assimilating the same data multiple times with increased measurement errors. MDA yields the same updated
mean and covariance as would be obtained from assimilating
the same data with the original measurement error covarianceand no iterations when errors are Gaussian, and all operators
are linear (Emerick and Reynolds 2012). For the nonlinear
case, EnKF-MDA partly resolves issues with nonlinearity andleads to smaller bias than a conventional EnKF. Emerick andReynolds (2013) developed the EnKS with MDA (EnKSMDA) for reservoir simulations, and Bocquet and Sakov(2014) showed IEnKS with MDA significantly outperforms
standard EnKF and EnKS in strongly nonlinear regimes with asimplified model. However, these methods have seen little attention in numerical weather prediction (NWP) communities.
While the convergence properties of these methods are unknown, numerical experiments performed by Evensen (2018)suggest they can provide accurate solutions for mildly nonlinear problems.
Last, particle filters (PFs) present a purely Bayesian filtering approach for state estimation, which avoids many ofthe linear/Gaussian assumptions of the above methods. PFs

VOLUME 149

provide a much more general, nonparametric estimate of
the model probability density function (PDF), which is advantageous for non-Gaussian problems as long as a sufficient number of ensemble members exist. Nevertheless,these methods can easily diverge when a relatively small
number of particles (ensemble members) are adopted fordata assimilation; see Bengtsson et al. (2008), Bickel et al.
(2008), and Snyder et al. (2008) for discussions on ensemble
size requirements for PFs. Several strategies are proposed toovercome this filter collapse and apply PFs to data assimilation problems for operational NWP models very recently.
One common effort to avoid filter divergence is to use localization, which restricts the influence of observations to
nearby state variables. For example, Poterjoy (2016) introduced the localized PF, which assimilates observations withindependent errors sequentially to combine sampled particles from a standard bootstrap PF with prior particles in amanner that satisfies a set of local constraints. Following thiswork, Poterjoy and Anderson (2016) and Poterjoy et al.
(2017, 2019) demonstrate that the local PF works well forhigh-dimensional systems. For these studies, the authors
compare the local PF with EnKFs for a simplified generalcirculation model and both idealized and real mesoscaleconvective systems in the Weather Research and Forecasting(WRF) Model, respectively. Even more recently, Potthast et al.
(2019) applied an alternative localized PF for global weather
prediction using the Icosahedral Nonhydrostic Weather andClimate (ICON) model, which marks the first successful test of a
PF in an operational framework. These studies provide an incentive to further explore the potential of localized PFs forweather prediction, especially considering the theoretical benefits they pose for assimilating remotely sensed measurements,
such as satellite radiance and radar reflectivity, which require
nonlinear measurement operators.
In addition to the methods described above, there are somenotable developments related to treatment of nonlinearity andnon-Gaussianity. For example, Bishop (2016) introduces the
GIGG-EnKF algorithm, which retains the accuracy of the
EnKF in the Gaussian case while lending it a high degree of
accuracy when the forecast and observation uncertainty are
gamma or inverse-gamma distributions. When conditions
are not suitable for EnKF, such as the distribution of theprior and observation are not Gaussian distribution, andthe observation operator is nonlinear, Amezcua and van
Leeuwen (2014) apply a preprocessing step known as
Gaussian anamorphosis to obtain state variables and observations that better fulfill the Gaussianity conditions.
Fletcher (2010) and Fletcher and Jones (2014) present variants of variational solvers for issues with lognormal andmixed lognormal Gaussian distributed background andobservation errors. While many methods have been proposed to deal with such difficult conditions, this study
mainly focuses on the tempered iteration approach, which isrelatively easy to implement in current NWP systems andcan deal with these problems well.
In this study, we discuss EnKF-MDA, EnKS-MDA, E4DVar,
4DEnVar, and the local PF data assimilation methods andtheir use in applications that require nonlinear measurement

Brought to you by University of Maryland, McKeldin Library | Unauthenticated | Downloaded 06/30/21 02:21 AM UTC

JULY 2021

2371

KUROSAWA AND POTERJOY

operators. We also examine the sensitivity of each method touser-specified parameters, which include ensemble size, covariance localization radius of influence (ROI), inflation coefficients, data assimilation window length (DAW), and thenumber of iterations and outer loops. The comparisons are
conducted with the 40-variable dynamical system introduced inLorenz (1996, hereafter L96), using numerical experiments
performed with conventional EnKF and EnKS techniques
as benchmarks. This study provides a necessary first stepin understanding the complexity of assimilating remotelysensed measurements in weather models, which will requireappropriate choices for data assimilation methodology going
forward.
Three main goals of these experiments are as follows: (1)investigate intrinsic limitations of current data assimilation
methodology for applications that require nonlinear measurement operators; (2) compare recently developed methods
designed for nonlinear/non-Gaussian applications with thosecurrently used for operational NWP; and (3) inform ongoing
efforts to design future geophysical modeling systems [e.g.,
NWP with Hurricane Analysis and Forecast System (HAFS)],
which will inevitably need to exploit remotely sensed
measurements.
The manuscript is organized in the following manner. In
section 2, we present algorithmic descriptions of each data
assimilation method. Section 3, describes settings for data assimilation experiments and results from the cycling experiments. Section 4 summarizes the main findings of this study anddiscusses the potential of the methods for real numerical
weather prediction.

xa 5 xf 1 K(y 2 Hxf ),(1)
21

K 5 Pf HT (HPf HT 1 R) ,T(2)P 5 (I 2 KH)P (I 2 KH) 1 KRK 5 (I 2 KH)P .
a

f

T

f

For the ensemble formulation, the covariance matrix P can be
defined as
P 5 EET ,(4)1
E 5 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ [dx(1) j . . . jdx(Ne ) ],Ne 2 1(5)where dx(l) is considered as a perturbation around x(l), which is
the lth member from an ensemble of Ne model states.
The Kalman filtering algorithm requires the computation of
Pa in (3). This process is equivalent to producing an appropriate analysis ensemble or ‘‘ensemble update,’’ which has asample covariance of Pa. For this study, all algorithms requiring an EnKF to update ensemble members use the serial ensemble square root filter (serial EnSRF; Whitaker and Hamill
2002). In general, this method provides a deterministic updateof the ensemble mean and perturbations about the ensemble
mean separately in a manner that satisfies the analysis meanand error covariance given by Kalman filter theory. The serial
EnSRF assumes an ensemble update of the following form:f~.
Ea 5 (I 2 KH)E

In this section, we present the mathematical framework foreach method, along with the dynamical system adopted for
performing numerical experiments. We use lowercase boldface
font to indicate vectors, uppercase boldface font to indicate
matrices, and italic font to indicate scalars and nonlinear
operators.
In this study, let x f be an N x -dimensional backgroundmodel forecast; let y be an N y -dimensional set of observations; let H be the tangent linear operator that converts
the model state to the observation space; let R be theN y 3 N y dimensional observation error covariance matrix; and let P be the N x 3 N x dimensional error covariance matrix. Superscript f and a denote forecast and analysis,
respectively.

a. EnKF
The EnKF is an approximate but efficient application of the
Kalman filter (Kalman 1960) and explicitly includes the time
evolution of error statistics, which operates effectively formoderately nonlinear dynamical systems. In EnKF, P is represented by ensemble members statistically. There is no needto consider the tangent linear model operator used in KF, so
EnKF has many advantages for nonlinear dynamics. The analyzed state xa is given by the following Kalman filter equations(e.g., Jazwinski 1970; Gelb et al. 1974):(6)

Andrews (1968) provides one solution, which involves Kalman
gain matrix for perturbations of the following form:~ 5 Pf HT [(HPf HT 1 R)21/2 ]T [(HPf HT 1 R)1/2 1 R1/2 ]21 .
K

2. Data assimilation methods

(3)

(7)

If observations are uncorrelated (R is diagonal), each observation is treated serially, which makes the terms HPfHT and R~ 5 aK
scalar. In this case, (3) can be simplified by assuming K
where a is a scalar value. The a was first derived by Potter(1964) assﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ !21
R
.
a5 11
f T
HP H 1 R(8)Thus, the serial version requires only the computation of ascalar factor to weight the traditional Kalman gain, andtherefore is no more computationally expensive than the
EnKF. In this study, observations are assumed to be independent of each other, which makes only the computation of(8) necessary. When assimilating a single observation throughthis formulation, K and H are vectors with Nx dimensions, andR is scalar. Therefore, for an individual observation, the terms
PfHT and HPfHT reduce to scalars and can be computed even if
the measurement operator is fully nonlinear, which is done by
applying this operator on each ensemble member before calculating sample statistics.

b. EnKSThe EnKS operates by storing ensemble members at pasttimes and then modifying them by a gain matrix that considers

Brought to you by University of Maryland, McKeldin Library | Unauthenticated | Downloaded 06/30/21 02:21 AM UTC

2372

MONTHLY WEATHER REVIEW

observations at the current time. Whitaker and Compo (2002)
introduced a serial ensemble square root smoother (serial
EnSRS), which uses Monte Carlo estimates of forecast–analysis
error cross covariances needed to compute the Kalmansmoother gain matrix. While they applied the serial EnSRS to
the fixed-lag Kalman smoother proposed by Cohn et al. (1994),
in this study, we apply it as a fixed-interval Kalman smoother.
Here, define a subscript notation mjn to indicate a quantity
at observation time m, which incorporates knowledge of allobservations up to and including time n. In this notation, (1)can be expressed as
xakjk 5 xfkjk21 1 K(y 2 Hxfkjk21 ) .

(9)

In the serial square root smoother, we use Pf(m,n) to denote across-covariance matrix between variables at times m and n.
The gain matrix K involves the forecast error cross-covariancematrix Pf(k,k2l) between xfkjk21 and xfk2ljk21 :K 5 Pf(k,k2l) HT (HPf HT 1 R)21 ,(10)

Pf 5 Efkjk21 EfkTjk21 ,

(11)

where

T
Pf(k,k2l) 5 Efkjk21 Efk2l
jk21 .

(12)

In the formulation of Cohn et al. (1994), this quantity is computed directly using the dynamical model because they developed
the fixed-lag smoother without ensembles. On the other hand, thefixed-lag smoother with ensembles uses the dynamical model only
when creating the background model forecast (Whitaker andCompo 2002). This idea can be directly implemented to the fixedinterval smoother. Note that the basic equations for the lag-0
implementation are identical to those of the serial EnSRF.

c. Multiple data assimilation (MDA)
Emerick and Reynolds (2012) introduced the MDA scheme,
which assimilates the same data multiple times using an inflated
covariance matrix of the measurement errors. They proved theequivalence between single and multiple data assimilations forthe linear-Gaussian case. Although MDA contains approximations for the fully nonlinear case and the equivalence does nothold for the nonlinear case, MDA benefits from the inclusion ofsmaller incremental ensemble corrections.
When the same set of observations are assimilated Na times,
the inflated measurement error covariance matrix is used in (2):K 5 Pf HT (HPf HT 1 ai R)21 ,(13)

the forecast and analysis state, the magnitude of the jump canbe overestimated by linear updates. This limitation of Gaussian
data assimilation techniques is observed frequently for the assimilation of all-sky radiance measurements in weather models,
which is one of the reasons to motivate the use of observation
error inflation (e.g., Minamide and Zhang 2017) and other ingenious approaches as described in section 1. By using an inflated
error covariance, a potentially large spurious update in the statevector is avoided. Going a step further, iterative techniques like
MDA replace single updates with a series of smaller updates,
which can correct filter or smoother updates that are too large.
In summary, the ensemble formulation of a fixed-interval
serial EnSRS, with and without MDA, are realized by thefollowing procedures. For DAW length l 5 0, the serial EnSRS
reduces to the serial EnSRF, and for Na 5 1, each iterative data
assimilation cycle with MDA reduces to a single-step data assimilation scheme, such as standard EnKF and EnKS.
Algorithm 1: EnKS with MDA cycle
1) Function MDAEnKS-MDA_cycle:
2) for t 5 1:time do
3)
if t is at the end of DAW then4)t0 ) t2l
5)
for i 5 1:iteration Na do6) for k 5 0:DAW length l do 7) xat0 jt0 1k ) Serial_EnSRS xft0 jt0 1k21 , xft0 1kjt0 1k21 , yt0 1k, ai RNaå 1/ai 5 1:

(14)

i51

Note that in this paper, we use ai 5 Na for i 5 1,. . .,Na for all
experiments with MDA. Rommelse (2009) and Emerick andReynolds (2012) suggest that when the assimilation of accurate
data in non-Gaussian regimes requires a ‘‘large jump’’ between

xft0 jt0 1k ) xat0 jt0 1kf
xt0 jt0 21 ) xat0 jt0 1l

8)
9)
10)
11)
12)
13)
14)

else

for m 5 1:Ne dof ( m)a(m)xt11jt ) Mxt0 jt

for i 5 1:iteration Na do
xatjt ) Serial_EnSRS(xftjt21 , xftjt21 , yt , ai R)fxtjt21 ) xatjt
for m 5 1: Ne do
f (m)
a(m)
xt11jt ) Mxtjt

15)
16)
17)
18) return

Algorithm 2: Serial EnSRS1) Function: Serial_EnSRS xft2kjt21 , xftjt21 , y, R2) for j 5 1: Ny do
h
i
1
f (1)
f (N )
3)
Eft2kjt21 5 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ dxt2kjt21 j . . . jdxt2kjet21
Ne 2 1hi
1f (1)f (N )
4)
Eftjt21 5 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ dxtjt21 j . . . jdxtjt21e
Ne 2 15)
Pf 5 Eftjt21 EftjTt21
6)
7)
8)

where

VOLUME 149

9)

Pf(t2k,t) 5 Eft2kjt21 EftjTt21
hi
K 5 Pf(t2k,t) H( j)T H(j) Pf H(j)T 1 R(j) 21

xat2kjt 5 xft2kjt21 1 K(y( j)2H( j) xftjt21 )0sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 121
R(j)
@
A
a5 11
H(j) Pf H(j)T 1 R(j)~ 5 aKK~ (j) Ef
Eat2kjt 5 Eft2kjt21 2 KH
tjt21
xat2kjt 5 xat2kjt 1 Eat2kjt
xft2kjt21 ) xat2kjt

10)
11)
12)
13)
14) return xat2kjt

Brought to you by University of Maryland, McKeldin Library | Unauthenticated | Downloaded 06/30/21 02:21 AM UTC

JULY 2021

2373

KUROSAWA AND POTERJOY

d. E4DVar and 4DEnVar
In this section, the equations of 4DVar, E4DVar, 4DEnVar
are introduced briefly. For further details on these methods, we
encourage readers to review the mathematical descriptions in
Liu et al. (2009), Poterjoy and Zhang (2015), and Bannister(2017). The 4DVar method seeks a solution that minimizes themisfit of a control variable to the background state xf0 at t 5 0 andobservations yt at times t 5 0, 1, 2,. . . , t. The minimization iscarried out with respect to increments dx0 from xf0 (Courtieret al. 1994). The cost function is expressed as the sum of background (Jb) and observation (Jo) terms:
J(dx0 ) 5 Jb (dx0 ) 1 Jo (dx0 )t

1
1
5 dxT0 B21 dx0 1 å (Ht Mt dx0 2 dt )T R21t (Ht Mt dx0 2 dt ) ,
2
2 t50
(15)
where B is the background error covariance and Mt is the tangent
linear model operator. The vector dt contains the innovations at
each time along a model trajectory from xf0 and is given by
dt 5 yt 2 Ht [Mt (xf0 )],(16)

where Mt and Ht are the nonlinear forecast model and observation operators, respectively. In practice, dx0 is replaced with
Uv, where v is the new control variable, and U is a square rootof the background error covariance matrix (B 5 UUT) (Lorenc
2003). The cost function in the control variable space and thegradient of the cost function with respect to the control variables become the following:
1
1 t
J(v) 5 vT v 1 å (Ht Mt Uv 2 dt )T R21t (Ht Mt Uv 2 dt ) ,
2
2 t50

(17)

t

=v J 5 v 1 å UT MTt HTt R21
t (Ht Mt Uv 2 dt ).

(18)

t50

For E4DVar and 4DEnVar, using a similar substitution
described above, dx0 is separated into two terms to include a
hybrid covariance in the variational cost function. For NWP
applications, the ensemble contribution of the hybrid covariance is often much greater than the static covariance (Kleist
and Ide 2015); however, such a choice is directly dependent
upon the quality of ensemble, ensemble size, and model error.
For the L96 model, Poterjoy and Zhang (2015) found the static
error covariance to have a major impact only when an imperfect model is used for data assimilation, which is not exploredin the current study. Therefore, we omit the use of a static error
covariance to reduce the number of parameters to examine forthis study. As a result, we havedx0 5 dxe0 5 Ue ve ,(19)

where dxe0 is the increment resulting from the ensembleestimated covariance. As described in Buehner (2005), Ue
can then be written as
Ue 5 [e(1) j . . . je(Ne ) ] ,(20)

P+C 5 Ue UeT ,sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
1
f (n)
(n)
e 5
3 diag(x0 2xf0 )C1/2 ,Ne 2 1(21)(n 5 1, 2, . . . , Ne ),(22)

where + indicates element wise multiplication, and C is thecorrelation matrix used for localizing the ensemble covariance.
From these equations, the cost function and the gradient of
E4DVar are found by substituting Ue for U and ve for v in (17)and (18). Using an ensemble forecast stored at each observation time in DAW, MtUe can be rewritten as
Mt Ue 5 [Mt e(1) j . . . jMt e(Ne ) ]
(N )

5 [^
e(t 1) j . . . j^et e ] ,sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
1
f ( n)^e(t n) 5
3 diag(xt 2 xft )C1/2
Ne 2 1sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
1
f (n)
5
3 diag(Mt (x0 )2 Mt (xf0 ))C1/2 .
Ne 2 1(23)(24)

By substituting (23) into (17) and (18), the 4DEnVar costfunction and the gradient can be expressed without the tangent
and adjoint model.
Note that while E4DVar uses tangent linear and adjoint
models to propagate a localized error covariance through the
DAW, 4DEnVar requires the localization of time covariances.
Most previous studies use the same correlation matrix at eachtime thus ignoring the complexity of introducing a localizationof time-dependent covariance (LTC) (Liu et al. 2009; Buehneret al. 2010; Liu and Xiao 2013; Fairbairn et al. 2014; Poterjoy
and Zhang 2015).
The method also allows for the use of either the nonlinear
operator Ht or the tangent linear operator Ht in its place.
This study explores both approaches in 4DEnVar experiments to identify which option presents the largest advantage for nonlinear operators. To perform the localization,
we calculate the tangent linear operator Ht at each time anduse it to propagate a localized error covariance through theDAW. Moreover, this study reruns the ensemble in outer
loops for 4DEnVar, despite the fact that it is prohibitively
costly for weather applications. This step is done to allow
for a more direct comparison with incremental E4DVar withouter loops.
To form a hybrid analysis, the variational solution is typically taken as the posterior mean and posterior perturbations
from an EnKF are recentered about this solution at the
middle of the time window (Zhang et al. 2009; Poterjoy et al.
2014). This approach is more consistent with the methodology adopted at major NWP modeling centers (Bannister2017). For the current study, we instead add posterior perturbations to the mean analysis at the end of each DAW. This
option has a number of advantages, namely, the EnKF assimilates measurements at the appropriate times over an assimilation window, thus providing an EnKF posterior mean that istheoretically equivalent to the 4DVar posterior mean in the

Brought to you by University of Maryland, McKeldin Library | Unauthenticated | Downloaded 06/30/21 02:21 AM UTC

2374

MONTHLY WEATHER REVIEW

absence of sampling error and nonlinearity. It also permits amore direct comparison of smoothers and filters explored inthis study.
In summary, the ensemble formulation of E4DVar and
4DEnVar are realized by the following procedures.
Algorithm 3: Ensemble/variational hybrid data assimilation
without static error covariance1) Function ensemble_variational_hybrid ðUe, xf0 , y, RÞ2) if 4DEnVar w/o LTC then3)Ht ) Ht4) while Outer Loop do5)
dt 5 yt 2 Ht [Mt (xf0 )]6)
while Inner Loop do
7)
switch Hybrid do
8)case E4DVar do
9)Dt ) MtUe
10)case 4DEnVar w/ LTC .or.
4DEnVar w/o LTC
h
i
(1)
(N )
et j . . . j^et e
11)Dt ) ^
1
12)
J(ve ) 5 veT ve
1 t 2e
1 åt50 (Ht Dt ve 2 dt )T R21t (Ht Dt v 2 dt )
2
t
T
T
21
13)
=ve J 5 ve 1 åt50 Dt Ht Rt (Ht Dt ve 2 dt )14)
ve 5 argmin(J(ve))15)
xf0 ) xf0 1 Ue vefa16) x0 ) x0
17) return xa0

e. The local PF
The current study uses the local PF proposed by Poterjoy
et al. (2019). For simplicity, this section highlights important
aspects of the local PF that are relevant to the comparisons performed in this study. Our experiments take advantage of additional regularization, tempering, and hybrid strategies that are
unique to the local PF, which are briefly discussed in this section.
For full details on this methodology, we refer readers to Poterjoy(2021, manuscript submitted to Quart. J. Roy. Meteor. Soc.).
The local PF assimilates observations serially, performing abootstrap PF update for particles projected onto the current
observation in the sequence, followed by a model-space update. For a given observation y, the model-space update replaces the standard bootstrap resampling step with one thatmerges sampled particles and prior particles:
xny 5 xy 1 r1 +(xkn 2 xy ) 1 r2 +(xkn 2 xy ),(25)

where xny is an updated particle, xn is the nth prior particle, xkn isthe nth sampled particle, xy is the localized posterior mean
based on importance weights that consider all observations up
to y, and r1 and r2 are derived to satisfy the posterior mean andvariance of marginals. The sampled particles are selectedfrom a bootstrap resampling of past updated particles using acumulative distribution formed by weights calculated fromparticle likelihoods for y. In general, the posterior particles
formed from linear combinations of the sampled and priorparticles are localized, because r1 and r2 are calculated basedon localized moments.

VOLUME 149

Poterjoy et al. (2019) provide several improvements to thePoterjoy (2016) local PF, which are aimed at preventing particle weight collapse. In addition, Poterjoy (2021, manuscript
submitted to Quart. J. Roy. Meteor. Soc.) introduces regularization and tempering methodology to further improve filter
performance when sampling error is large. In short, regularization raises particle weights to a power b, which is predetermined to yield marginal particle weights that have a
specified ‘‘effective sample size,’’ similar to the methodology described in Poterjoy et al. (2019). Regularization acts
as a heuristic means of preventing weight collapse, similar toobservation error inflation. It provides a strategy for assimilating observations through tempered iterations (Neal
1996), each with a unique set of b coefficients. Unlike regularization, tempering does not introduce bias in the posterior estimate.
The method also benefits from the use of a mixing parameter, g, to increase particle diversity in the vicinity of observations. As described in Poterjoy (2021, manuscript submitted
to Quart. J. Roy. Meteor. Soc.), r1 in (25) is multiplied by g,which introduces a smooth ‘‘jittering’’ of particles. The coefficients in r2 are then modified so that the first two posterior
moments are still maintained.

3. Cycling data assimilation experiments
We perform separate sets of data assimilation experiments
to investigate limitations for nonlinear applications and examine the sensitivity of the methods to user-specified parameters. These parameters include the number of iterations,
DAW, ensemble size, ROI, inflation, and measurement operators. The first two sets of experiments focus primarily on key
parameters for smoothers, which are known to be sensitive tononlinearity in model dynamics and measurement operators.
These parameters are the number of iterations and DAWlength. The third set of experiments focuses more broadly on thecomparison between filters and smoothers. For this purpose, we
select three types of observation networks, each differing primarily in choice of measurement operator. The system parameters for each of these cases are summarized in Table 1.

a. Experimental design1) MODEL
We examine several aspects of the data assimilation methods
by performing idealized numerical experiments with the L96
model (Lorenz 1996; Lorenz and Emanuel 1998). The model
consists of variables xi for i 5 1, 2,. . . , Nx, which are equally
spaced on a periodic domain. The variables are evolved in timeusing the set of differential equations:dxi
5 (xi11 2 xi22 )xi21 2 xi 1 F ,
dt

(26)

with cyclic boundaries: xi1Nx 5 xi and xi2Nx 5 xi . We integrate(26) forward numerically using the fourth-order Runge–Kutta
method with a time step of 0.05 [units defined arbitrarily as 6 h;
see Lorenz (1996)]. For this study, we fix Nx at 40 and use
F 5 8.0, which causes the model to behave chaotically.

Brought to you by University of Maryland, McKeldin Library | Unauthenticated | Downloaded 06/30/21 02:21 AM UTC

JULY 2021

2375

KUROSAWA AND POTERJOYTABLE 1. Configuration of cycling data assimilation experiments.

Expt

H(x)

sy

Ne

ROI

a,g

Ny

Dt (h)

Linear Case
Nonlinear Case 1
Nonlinear Case 2xx2
log(jxj)

1.0
1.0
0.1

10, 4010, 4040, 100

Variable
Variable
Variable

Variable
Variable
Variable20
20, 15, 10
20

6
6, 24
6

2) OBSERVATIONSIn this study, we create observation networks of Ny 5 10, Ny 5
15, and Ny 5 20 observations that are evenly spaced on model
grid points. Note that for the case Ny 5 15, we line up the observation points so that they were evenly distributed (i.e., at grid
points 1, 4, 6, 9, 12, 14, 17, 20, 22, 25, 28, 30, 33, 36, 39). We
simulate measurements every time step (6 h) by selecting values
from a truth simulation, applying one of the operators discussed
below, then adding uncorrelated Gaussian errors selected fromN(0, s2y I), where s2 is the measurement error variance.
Experiments include three forms of measurement operator.
The ‘‘Linear Case’’ uses an H that selects model variables to be
directly observed; i.e., H(x) 5 ^x, where ^x is a subset Ny variables in x chosen by H. The ‘‘Nonlinear Case 1’’ extends H tobe quadratic: H(x) 5 ^x+^
x. The ‘‘Nonlinear Case 2’’ introduces
log and absolute value operators to the interpolated values:
H(x) 5 log[ABS(^x)], where ABS indicates the absolute valueof each element. The second and third operators produce weak
and strong nonlinearities, respectively. Note that we apply asimple gross error check for the third measurement operator toprevent observations from being assimilated if the value ofABS(^x) is extremely small. If the difference between the observed value and the background state is greater than four times
the standard deviation of the background state, the observation
will not be assimilated. Observation error standard deviations
are set to sy 5 1.0 for the first two experiments, but reduced tosy 5 0.1 for the third case to compensate for the smaller information content provided by this observation network.

3) OBSERVATION TIME LINE AND VERIFICATION
Observations are assimilated over a 3650-day period, androot-mean-square errors (RMSEs) from the last 3550 days areused to quantify the accuracy of the posterior analyses. The
first 100 days of data assimilation act as a spinup period to allow
members time to reach quasi-steady posterior solutions for the
given setup of the model and observation network.
In the first sets of experiments described below, we perform
direct comparisons of the different smoothers used for thisstudy. For these experiments, we calculate RMSEs at the beginning of the DAW (smoother solution), because it more directly indicates how much information is being extracted fromobservations at future times. For experiments shown later in thissection, which compare different forms of smoothers and filters,
we calculate RMSEs at the end of the DAW (filter solution).

4) TREATMENT OF SAMPLING ERRORS
Potential sources of bias in the estimation of the posterior
include small ensemble sizes relative to the state dimensions,model errors, nonlinearities, and assumptions used to form
data assimilation algorithms. Therefore, heuristic covariancelocalization strategies are needed to reduce noise introduced
from ensemble error approximations by performing a Schurproduct between this matrix and an empirically defined correlation matrix with a tunable length scale parameter, or ROI.
For this purpose, we use the fifth-order correlation function
given by Eq. (4.10) of Gaspari and Cohn (1999).
The posterior covariance is inflated by replacing ensemble
perturbations with linear combinations of posterior and priorperturbations, which is known as a covariance relaxation
method (Zhang et al. 2004):
0a

0a

0

xn ) (1 2 a)xn 1 axnf .

(27)

The a in (27) is called the ‘‘relaxation coefficient’’ and ranges
from 0 to 1, where a 5 0 implies no inflation. We adopt thisinflation strategy to remain consistent with Poterjoy and Zhang(2015), who perform a similar comparison of ensemble data
assimilation algorithms, including hybrid covariance forms of
E4DVar and 4DEnVar.
As previously stated, the local PF uses a mixing parameter tomaintain particle diversity during updates. While this approachis effective at preventing filter divergence with small ensembles, it does not directly increase prior or posterior error variance in the same manner as relaxation. Similar to the a used in
the relaxation method the coefficient g is a scalar between 0
and 1. It further mixes prior particles and resampled particles
everywhere particles are updated in state space, including in
the vicinity of measurements.

b. Results
1) SENSITIVITY TO THE NUMBER OF OUTER ITERATIONS
The variational and MDA techniques present different iterative strategies for coping with nonlinearity in model dynamics and measurement operators. For the first set ofexperiments, we explore the sensitivity of these methods to thenumber of iterations. In addition to providing a direct comparison of different smoothers for a nonlinear application,these experiments help motivate choices for iteration numberin the filter/smoother comparisons that follow. As previously
stated, we also explore the advantage of LTC, which is a localization of the ensemble covariance at each observation timein the window calculated with the tangent linear operator Ht ateach time for nonlinear operators.
Figure 1 shows mean RMSEs of EnKS-MDA, E4DVar,
4DEnVar with LTC, and 4DEnVar without LTC from experiments with Nonlinear Case 1. Ensemble size Ne, relaxation
coefficient a, and DAW are fixed at 10, 0.3, and 24 h, respectively. We find this window length to be sufficient for exploring
sensitivity to outer loops without adding computational cost.
We do not show results using Nonlinear Case 2 because all

Brought to you by University of Maryland, McKeldin Library | Unauthenticated | Downloaded 06/30/21 02:21 AM UTC

2376

MONTHLY WEATHER REVIEWVOLUME 149

FIG. 1. Mean analysis RMSEs as a function of the number of iteration or outer loop. Results
are shown for the Nonlinear Case1. Values are from the experiment with EnKS-MDA (triangles), E4DVar (circles), 4DEnVar without LTC (squares), and 4DEnVar with LTC (diamonds), and ROI set to 1 (blue), 3 (red), and 5 (green). The RMSEs are calculated at the startof the DAW (smoother solution).

methods tested in this study (other than the PF) experiencefilter divergence when measurements are simulated with thisoperator. These results are discussed in the filter/smoothercomparisons below.
For the observation networks tested in this study, we findthat increasing the number of iterations has little impact onmean error for EnKS-MDA. For E4DVar and 4DEnVar,
however, we confirm that multiple outer loops are required for
optimal performance. Under various circumstances, outer
loops are also needed to prevent filter divergence with thenonlinear measurement operator. For example, E4DVar with
ROI fixed at 1 and a single outer loop shows a worse score thanwith multiple iterations. We also find that the minimum number of outer loops required to prevent filter divergence issensitive to ROI. E4DVar experiments using an ROI of 3 and 5
require 2 and 3 outer iterations, respectively. Nevertheless, theimprovements of multiple iterations beyond these numbers
become negligible once a sufficient number is reached.
We also find E4DVar to be more stable than 4DEnVar forthe tested observation networks. Recall, this method uses thetangent linear model to propagate increments along a nonlinear trajectory to future times, and its adjoint to propagate
sensitivity gradients backward from observation times to thebeginning of the DAW. The trajectory is updated between
outer iterations to ensure that values propagated by thetangent linear and adjoint remain small enough for linear
approximations to remain valid. In addition, the input of
ensemble error covariance at a single time in this process (atthe beginning of the DAW) greatly simplifies the removal of
spurious error correlations through localization (Fairbairn
et al. 2014; Poterjoy and Zhang 2015). For this reason, we findconfigurations of 4DEnVar that use LTC to be more stable
than configurations without LTC. Based on this finding, we
use this strategy for all remaining 4DEnVar experiments.

2) SMOOTHER PERFORMANCE AS A FUNCTION OF DATAASSIMILATION WINDOW LENGTH

Several of the methods examined in this study are smoothers,
which are sensitive to the choice of DAW. For the next set ofexperiments, we compare mean RMSEs of EnKS, EnKS-MDA,
E4DVar, and 4DEnVar as a function of DAW (Fig. 2). As statedabove, the verification for these experiments focuses on theposterior smoothing density; i.e., the analysis at the beginning of
the DAW. For these experiments, we fix the ensemble size Ne,
relaxation coefficient a, and ROI at 10, 0.3, and 3, respectively.
The number of iterations (MDA) and outer loops (Var) are both
set to 3. These decisions are based on results from the previous
set of experiments, showing little benefit beyond 3 iterations forchosen model and observation networks. As we revisit later, inexperiments with the Nonlinear Case 1, the observation value iscloser to the truth all the time, making the RMSEs and order ofmagnitude smaller than in the Linear Case.
We start by examining the impact of MDA on the EnKS.
Our experiments show that MDA provides slight benefits
over noniterative configurations, even at DAW length l 5 0 h
and linear H (Fig. 2a). Note that EnKS is identical EnKF forthis DAW length, so no benefits are expected from the iterations. One possible reason for the difference in skill between
EnKS and EnKS-MDA at DAW length l 5 0 h is due to small
differences in how ensemble perturbations are adjusted
through iterative steps. For linear cases with Gaussian prior,
MDA yields the same posterior mean and covariance as
would be obtained without iterations. As suggested by
Rommelse (2009), the extra uncertainty included in measurements during each iteration ensures that adjustmentsfrom prior to posterior values are dampened, which is beneficial when linear updates overestimate the true impact ofmeasurements that relate nonlinearly to model variables.
Therefore, MDA provides an opportunity for the EnKF toremove overadjustments that may occur during previous iterations. We suspect that a combination of serial processingof observations and iterative updates of members leads toslight improvements in how the EnKF samples from theposterior density, which is assumed to be non-Gaussian because of the nonlinear model. This finding explains why the
MDA approach yields small improvements in posterior estimates over successive data assimilation steps, which is also
explored later.

Brought to you by University of Maryland, McKeldin Library | Unauthenticated | Downloaded 06/30/21 02:21 AM UTC

JULY 2021

KUROSAWA AND POTERJOY

2377

FIG. 2. Mean analysis RMSEs as a function of smoother lag. Results are shown for (a) Linear
Case and (b) Nonlinear Case 1, with the EnKS (blue), the EnKS-MDA (red), the E4DVar(green), and the 4DEnVar (magenta). The number of iterations and outer loops is fixed at 3 forboth cases. The RMSEs are calculated at the start of the DAW (smoother solution).

The advantage of the EnKS-MDA over the EnKS with the
DAW length l . 0 h is shown in both the linear and nonlinear
cases. For both experiments, the MDA scheme resolves issues
with the nonlinearity of the model and observation measurement operators in DAW. EnKS is stable even with thelonger DAW, but the quality of the analysis starts to degrade as the DAW length is increased beyond a certain
point, because sampling error increases as the DAW become longer. Compared to 4DEnVar, EnKS is more stablewith longer DAW. This indicates that the forecast error
covariance matrix used for smoother is approximated moreaccurately by cross-covariance matrix (Pf(k,k2l) ) in EnKS
than by ensemble-based error covariance in 4DEnVar.
Unlike the variational methods, the EnKS samples directlyfrom the smoothing density rather than using a hybrid
strategy of recentering EnKF perturbations about a variational solution. Furthermore, the 4DEnVar experimentcontains higher RMSEs than E4DVar because of the difficulty required in removing sampling errors from temporalerror covariances when Ne is small (Fairbairn et al. 2014;
Poterjoy and Zhang 2015).

3) FILTER PERFORMANCEIn this section, we present results from experiments thatexamine the sensitivity and limitations of EnKF, EnKF-MDA,
EnKS, EnKS-MDA, E4DVar, 4DEnVar, and the local PF to
ROI, relaxation coefficient a, PF mixing coefficient g, and theobservation measurement operators. For all experiments,
DAW for EnKS, EnKS-MDA, E4DVar, and 4DEnVar is set
to 24 h, and the number of iterations and outer loops are set to3. For the local PF, the regularization operates only when theeffective ensemble size Neff falls below a target value of Neff.
t
5 0:5 3 Ne for all experiments.
The target Neff is fixed at Neff
We define filter divergence objectively by flagging configurations that produced 100-cycle average RMSEs larger than 2with NA for ‘‘not available’’ in the figures.
Figure 3 shows mean RMSEs from the experiment with theLinear Case. Results from all methods, which use a fixed ensemble size Ne of 10, are displayed in charts that show RMSE
as a function of tunable variables used to reduce the impact ofsampling error. For example, Fig. 3 demonstrates that the
optimal ROI and a are comparable for EnKF, EnKF-MDA,
EnKS, EnKS-MDA, E4DVar, and 4DEnVar. In most cases,
the optimal scores are typically found near values that lead tofilter divergence. RMSEs from the local PF are slightly worse
due to the small number of particles used in these experiments.
Figure 4 shows results from experiments with the same settings
except Ne is increased to 40. As expected, all methods become
more stable and require less localization (larger ROI) and lessinflation (smaller a and g) as Ne increases. Comparing the results of the local PF from Figs. 3 and 4, it is clear that the larger
ensemble size is required for the local PF to outperform themethods with a Gaussian prior with the tested observation
network. EnKS shows clearly better performances than EnKF,

Brought to you by University of Maryland, McKeldin Library | Unauthenticated | Downloaded 06/30/21 02:21 AM UTC

2378

MONTHLY WEATHER REVIEWVOLUME 149

FIG. 3. Mean analysis RMSEs estimated for a range of (a)–(f) relaxation coefficient a and (g) PF mixing coefficient g and ROI. Results are shown for experiments with the Linear Case and ensemble size is fixed at 10. Blackshading indicates higher RMSEs, NA indicates that filter divergence occurs during the experiment, and the smallest
errors are indicated by the black boxes. The RMSEs are calculated at the end of the DAW (filter solution).

and MDA makes EnKF and EnKS slightly improved, evenwith a linear measurement operator because of the reason
mentioned in section 3b(3).
Results from Nonlinear Case 1 experiments using Ne 5 10
are shown in Fig. 5. Unlike experiments with the Linear operator, filter divergence occurs without setting strict limits on
ROI and inflation coefficients for all methods. Despite thenonlinear measurement operator in these experiments, we findno benefits from the assimilation methods designed specificallyfor non-Gaussian applications, namely, EnKF-MDA and the
local PF. We believe this result occurs because of the accuracy
and frequency at which these measurements are collected. For
model variables that can reach magnitudes of O(10), measuring the square of these variables with an error variance of 1

Brought to you by University of Maryland, McKeldin Library | Unauthenticated | Downloaded 06/30/21 02:21 AM UTC

JULY 2021

KUROSAWA AND POTERJOY

2379

FIG. 4. As in Fig. 3, but for ensemble size fixed at 40.

yields highly accurate information for characterizing the posterior. This factor, combined with the frequency of thesemeasurements lead to prior and posterior members that remain close to the truth at all times, thus making Gaussian assumptions more valid. We revisit this property of the Nonlinear
Case 1 measurement operator in the next section.
These experiments also continue to show clear benefits of
E4DVar and 4DEnVar over EnKF, both in terms of stability
and accuracy. We hypothesize that the 4D data assimilation
methods are less sensitive to sampling noise, which becomes

the dominant source of bias in mildly nonlinear regimes.
Likewise, we find E4DVar to be more stable than 4DEnVar
when Ne is small, owing mostly to the localization strategy
adopted by this method. We note that all algorithms approachsimilar RMSEs as ensemble size increases; i.e., Fig. 6 shows
results with Ne 5 40 for the same observation network. The
reason why E4DVar and 4DEnVar are more stable than EnKS
is due to the small number of ensembles and the nonlinear
observations that prevent from accurately estimating of thecross-covariance matrix in the Serial EnSRS.

Brought to you by University of Maryland, McKeldin Library | Unauthenticated | Downloaded 06/30/21 02:21 AM UTC

2380

MONTHLY WEATHER REVIEWVOLUME 149

FIG. 5. Mean analysis RMSEs estimated for a range of (a)–(f) relaxation coefficient a and (g) PF mixing
coefficient g and ROI. Results are shown for experiments with the Nonlinear Case 1 and ensemble size is fixedat 10. Black shading indicates higher RMSEs, NA indicates that filter divergence occurs during the experiment, and the smallest errors are indicated by the black boxes. The RMSEs are calculated at the end of theDAW (filter solution).

Figure 7a shows the mean RMSEs from experiments of thelocal PF that use measurements simulated with Nonlinear Case2 and Ne 5 40. For this configuration, filter divergence occurs inall methods except the local PF, owing to the strong nonlinearity in the measurement operator. This observation network
presents a case where nonlinearity in the application becomes amuch larger factor than sampling error in ensemble-estimatedprior and posterior distributions. Even with Ne 5 100, theGaussian-based methods fail to provide stable solutions despite the potentially large amount of information contained inthese measurements, as indicated by the low RMSEs in thelocal PF posterior (Fig. 7b). Since the local PF makes no

Brought to you by University of Maryland, McKeldin Library | Unauthenticated | Downloaded 06/30/21 02:21 AM UTC

JULY 2021

KUROSAWA AND POTERJOY

2381

FIG. 6. As in Fig. 5, but for ensemble size fixed at 40.

parametric assumptions about prior densities, non-Gaussian
observation-space priors, which are produced by nonlinear
measurement operators, do not have a negative impact on thefilter. Therefore, it can continue to extract information fromthe observation network regardless of nonlinearity in H.
These results confirm past studies, showing that local PF
provides benefits when Ne is sufficiently large or when theobservation operator is strongly nonlinear. It also demonstrates limitations in iterative techniques for cases where theobservation function is quadratic and the posterior may bebimodal.

4) FILTER PERFORMANCE FOR SPARSE OBSERVATIONNETWORKS

Using the mildly nonlinear observation operator (Nonlinear
Case 1), we investigate the behavior of each method for increasingly sparse observation networks. These experiments use
an observation frequency of 24 h, which is increased from 6 h in

Brought to you by University of Maryland, McKeldin Library | Unauthenticated | Downloaded 06/30/21 02:21 AM UTC

2382

MONTHLY WEATHER REVIEWVOLUME 149

FIG. 7. Mean analysis RMSEs estimated for a range of PF mixing coefficient g and ROI. Filter divergence occurs
in all methods except the local PF, so only results of the local PF are shown for experiments with the Nonlinear Case2 and ensemble size is fixed at (a) 40 and (b) 100. Black shading indicates higher RMSEs, NA indicates that filter
divergence occurs during the experiment, and the smallest errors are indicated by the black boxes. The RMSEs arecalculated at the end of the DAW (filter solution).

previous experiments, and Ny 5 20, 15, and 10 for equally
spaced measurements at each observation time. We also fix theDAW for smoothers at 48 h; see Table 1 for full summary.
These results are summarized in Figs. 8–10 using the samegraphics adopted in the previous section comparing filter
performance.
Compared to EnKF, the performance of EnKS becomesslightly worse for these observation networks. As discussed in
Evensen and van Leeuwen (2000), the EnKS differs from the
EnKF by computing updates of the model parameters using all
the observations in DAW simultaneously rather than using
recursive updates in time. Therefore, with these settings, therecursive updates of EnKF keep the model solutions close to
the truth at any given time during the experiment, and operate
on marginal densities that are relatively close to Gaussian at
any given time. While posterior marginals of the smoothingdensity are expected to be close to Gaussian at the beginning of
the DAW (Morzfeld and Hodyss 2019), marginals near the endof the DAW can evolve non-Gaussian characteristics becauseof nonlinearity in the model.
The benefits of MDA for EnKF are clearly shown in Figs. 8
and 9. For suboptimal configurations of the EnKF, priormembers exhibit a larger variance thus allowing nonlinearity in
H to become a significant source of bias for Gaussian methods.
Therefore, the optimal EnKF configuration remains almost thesame with MDA, but the set of parameters over which the filter
remains stable becomes larger than that of the standard EnKF.

For these observation networks, careful choices of ROI anda are sufficient for mitigating bias caused by Gaussian assumptions, but MDA helps prevent filter divergence whenthese parameters are improperly chosen.
For a long DAW (48 h) E4DVar becomes more stable than
4DEnVar with Ny 5 20 (Fig. 8), but both methods diverge
when observation density is decreased further (Figs. 9 and 10).
For these experiments, we find EnKS-MDA to be more accurate than the EnKS and much more stable than the variational methods. This result is anticipated in nonlinear regimes,
since incremental updates reduce potential overadjustments by
the ensemble smoother over the time window. As previously
stated, the improved performance over E4DVar and 4DEnVar
for sparse observation networks (Fig. 9) must follow from theability of EnKS-MDA to sample directly from the smoothingdensity, rather than relying on a hybrid approach, which is aclear advantage of this method. Algorithmically, the EnKS
operates in a manner that is very similar to 4DEnVar, but withthe added benefit of updating ensemble perturbations about
the posterior mean, rather than recentering EnKF perturbations about the posterior mode.
For the experiment with Ny 5 15, we also verify the second
moment of the posterior to examine potential shortcomings in
uncertainty estimates. The observation network and ensemble
size used in these simulations poses challenges for several data
assimilation method used here, in that filter divergence isprevented for a narrower range of parameters than previous

Brought to you by University of Maryland, McKeldin Library | Unauthenticated | Downloaded 06/30/21 02:21 AM UTC

JULY 2021

KUROSAWA AND POTERJOY

2383

FIG. 8. As in Fig. 6, but for the frequency of observations and DAW fixed at 24 and 48 h, respectively.

experiments. Figure 11 shows the ratio of spread to RMSEs,
indicating whether the ensemble spread is overestimated orunderestimated with respect to the RMSE. The results of all
methods are presented except E4DVar and 4DEnVar, which
do not estimate posterior variance—recall that ensemble perturbations are updated using an EnKF instead. Ideally, thespread and RMSE should be equivalent, but sampling error
and assumptions made during data assimilation may lead toinconsistent results. Likewise, heuristic techniques for treating

sampling errors, such as localization and covariance relaxation
can also introduce suboptimal uncertainty estimates. For all
filters and smoothers examined in this study, the best match
between spread and RMSE tends to occur when RMSE is at aminima (Figs. 11a–d). The farther away from the optimal parameter settings, the larger the mismatch between spread and
RMSE. As such, filter divergence occurs when the spread begins to become overestimated or underestimated for all
methods (Figs. 11a–d).

Brought to you by University of Maryland, McKeldin Library | Unauthenticated | Downloaded 06/30/21 02:21 AM UTC

2384

MONTHLY WEATHER REVIEWVOLUME 149

FIG. 9. As in Fig. 8, but for the number of observations fixed at 15.

Despite the difficulty posed by these observation networks, we find that the local PF can be configured to produce
stable results, even for data-sparse regimes, which was expected for this method (Poterjoy 2021, manuscript submitted to Quart. J. Roy. Meteor. Soc.). This property of the local
PF is illustrated for the Ny 5 10 case, where it is the only
method that does not diverge for all parameter value(Fig. 10). These results demonstrate challenges that exist forthe mildly nonlinear observation operator as the spatial andtemporal density of measurements decreases to yield larger
prior uncertainty.

5) LOCAL PF PERFORMANCE AS A FUNCTION OFENSEMBLE SIZE

Figure 12 shows the mean RMSEs of the local PF as afunction of ensemble size. These experiments use a fixed PFtvalues of 0.2 3 Ne andmixing coefficient of g 5 0.3 and two Neff
0.8 3 Ne. The results are similar for the cases with the linear

Brought to you by University of Maryland, McKeldin Library | Unauthenticated | Downloaded 06/30/21 02:21 AM UTC

JULY 2021

KUROSAWA AND POTERJOY

2385

FIG. 10. As in Fig. 8, but for the number of observations fixed at 10. Filter divergence occurs in all methods except
the local PF, so only results of the local PF are shown.

and mildly nonlinear measurement operators (Figs. 12a,b) inthat optimal ROI increases with ensemble size.
This is because the large ensemble size yields fewer sampling
errors thus needing less localization. For the strongly nonlinear
measurement operator, however, the difference in RMSEs forthe range of ROI choices become small as the ensemble size
increases (Fig. 12c). This result may reflect either the limited
information contained in these measurements. That is, because
they only observe the log of the absolute value of variables,
distant multivariate updates from these measurements are
truly very small, thus requiring very large ensemble sizes toestimate accurately. They may also suggests that sampling errors and other factors, such as assumptions made by local PF
update equations, become less dominant for nonlinear applications of this type.
Furthermore, the experiments demonstrate a dependence ofton ensemble size. When the ensemble size is small,
optimal Nefftshow more accurate results. As
experiments with higher Nefftthe ensemble size increases, the lower Neffshows smallert
can result
posterior RMSEs. This suggests that the larger Neffin overinflation when the ensemble size is large.

4. Conclusions
In geophysical models, such as those used for numerical
weather prediction, strongly nonlinear model dynamics andmeasurement operators can cause data assimilation methods tobe biased. This study examines several procedures that aredeveloped to overcome challenges posed by nonlinear operators, such as periodic relinearization of tangent linear and adjoints in variational schemes, likelihood factorizations adopted
by iterative ensemble filters and smoothers, and localized
particle filters. These methods—some of which were originallydesigned for applications outside the weather community—arecompared with methods currently used for operational NWP,
namely, EnKFs and hybrid variational methods with andwithout model adjoints.
This study adopts the 40-variable model of Lorenz (1996) toexamine the selected data assimilation approaches. The small
dimension of this model allow for extensive testing of eachtechnique using a large variety of observation networks, eachvarying in density and the type of observations provided. For

several observation networks used in this study, relinearizationof the model and measurement operators between outer iterations are required to prevent filter divergence. Once a sufficient number of outer iterations are reached to achieve stable
results, the improvements are negligible.
The wide range of observation networks examined in thisstudy yields a diverse set of results, which are summarizedusing posterior RMSEs. We acknowledge that this metric is not
ideal for non-Gaussian regimes, particularly those characterized as multimodal. Nevertheless, the sharp failure of various
techniques for non-Gaussian problems are easily identified bylarge values of RMSEs.
Each method examined in this study has clear advantages for
specific regimes—which are identified to be a function ofsampling error, nonlinearity in measurement operators, andobservation density. This finding motivates the use of different choices of data assimilation methodology, depending
on application.
The ensemble-variational smoother with an adjoint model,
E4DVar, produces smaller RMSEs than 4DEnVar for all observation networks tested in this study. It also outperforms allother methods in regimes where sampling error is high, but themodel solution is well-constrained by numerous accurate
measurements; i.e., in weakly nonlinear regimes. This study
also compares variational methods to an ensemble smoother,which is adapted from the fixed-lag EnSRS of Whitaker andCompo (2002). For regimes where sampling error is a more
dominant sources of posterior bias than nonlinearity, the EnKS
performs better than its filter counterpart. Adding iterations to
EnKF and EnKS updates through MDA results in improved
results for all nonlinear regimes, particularly for sparse observation networks and long DAW lengths. The EnKS with
MDA is also found to outperform all methods for data assimilation problems characterized by high sampling error andweak nonlinearity. Likewise, it provides stable results in nonlinear regimes that cause E4DVar and 4DEnVar to experiencefilter divergence. For applications of this type, EnKS-MDA
benefits from its ability to sample directly from the posteriorsmoothing density, rather than relying on a separate EnKF to
update perturbations about a maximum likelihood solution.
Furthermore, ensemble filters outperform smoothers whennonlinearity in measurement operators or model dynamics

Brought to you by University of Maryland, McKeldin Library | Unauthenticated | Downloaded 06/30/21 02:21 AM UTC

2386

MONTHLY WEATHER REVIEWVOLUME 149

FIG. 11. Ratio of ensemble spread to mean analysis RMSEs estimated for a range of (a)–(d) relaxation coefficienta and (e) PF mixing coefficient g and ROI. The experimental setting is as in Fig. 9. NA indicates that filter divergence occurs during the experiment. The RMSEs and spread are calculated at the end of the DAW (filter
solution).

have a dominant role in the data assimilation applications.
This finding is consistent with past studies that compare filtersand smoothers for problems of this type (Evensen and van
Leeuwen 2000). For highly nonlinear regimes, the local PF is
the only method that produces accurate results. The benefit of
PF-based methodology, however, comes with the trade-off ofbeing more sensitive to sampling error. Therefore, it requires
large ensemble sizes to produce RMSEs as low as ensemble
and variational smoothers for quasi-linear regimes.
Owing to the nature of this study, all comparisons are performed in an idealized framework. These findings will ultimately help guide future data assimilation decisions for real
geophysical problems, where the computational cost of exploring the sensitivity of data assimilation methodology andparameters is prohibitive. The major findings of this study

demonstrate when to expect Gaussian filters and smoothers tobe suboptimal and under what conditions iterative techniques
provide added value over conventional methods. Choices ofnonlinear measurement operators in this study are motivatedby challenges faced by high-impact weather events, such as
severe convective storms and tropical cyclones. In particular, all-sky satellite radiance measurements provide extensive, near-continuous data coverage for tropical cyclones
over open oceans. These measurements are often difficult touse, owing to the highly non-Gaussian (often multimodal)
observation-space priors produced by nonlinear measurement operators. New operational weather prediction systems, such as NOAA’s Hurricane Analysis and Forecast
System (HAFS), will ultimately need to overcome barriers
that currently exist in Gaussian-based data assimilation

Brought to you by University of Maryland, McKeldin Library | Unauthenticated | Downloaded 06/30/21 02:21 AM UTC

JULY 2021

KUROSAWA AND POTERJOY

2387

FIG. 12. Mean analysis RMSEs of the local PF as a function of ensemble size. Results areshown for (a) Linear Case, (b) Nonlinear Case 1, and (c) Nonlinear Case 2. Values are from theexperiment with Neff fixed at 0.20 3 Ne (solid lines) and 0.80 3 Ne (dashed lines), and ROI fixedat 2 (blue), 5 (red), and 8 (green).

methodology to fully leverage measurements of this type, as
several operational centers have made significant advancements
to cope with the difficult conditions in the past years. Experiments
performed in this study motivate applications of iterative ensemble approaches and the local PF for problems of this type.
Acknowledgments. Funding for this work was provided by
NOAA Grant NA20OAR4600281 and NSF/CAREER Award
AGS1848363.
Data availability statement. All software used to generate
results for this study is available upon request from thecorresponding author.
REFERENCES
Amezcua, J., and P. J. van Leeuwen, 2014: Gaussian anamorphosis
in the analysis step of the EnKF: A joint state-variable/
observation approach. Tellus, 66A, 23493, https://doi.org/
10.3402/tellusa.v66.23493.
Andrews, A., 1968: A square root formulation of the Kalman covariance equations. AIAA J., 6, 1165–1166, https://doi.org/
10.2514/3.4696.
Bannister, R. N., 2017: A review of operational methods of variational and ensemble-variational data assimilation. Quart.
J. Roy. Meteor., 143, 607–633, https://doi.org/10.1002/qj.2982.
Bengtsson, T., P. Bickel, and B. Li, 2008: Curse-of-dimensionality
revisited: Collapse of the particle filter in very large scale
systems. Probability and Statistics: Essays in Honor of David
A. Freedman, D. Nolan and T. Speed, Eds., Vol. 2, Institute ofMathematical Statistics, 316–334.
Bickel, P., B. Li, and T. Bengtsson, 2008: Sharp failure rates for thebootstrap particle filter in high dimensions. Pushing the Limitsof Contemporary Statistics: Contributions in Honor of Jayanta
K. Ghosh, B. Clarke and S. Ghosal, Eds., Vol. 3, Institute ofMathematical Statistics, 318–329.

Bishop, C. H., 2016: The GIGG-EnKF: Ensemble Kalman filteringfor highly skewed non-negative uncertainty distributions.
Quart. J. Roy. Meteor. Soc., 142, 1395–1412, https://doi.org/
10.1002/qj.2742.
Bocquet, M., and P. Sakov, 2014: An iterative ensemble Kalmansmoother. Quart. J. Roy. Meteor. Soc., 140, 1521–1535, https://doi.org/10.1002/qj.2236.
Bonavita, M., P. Lean, and E. Hólm, 2018: Nonlinear effects in
4D-Var. Nonlinear Processes Geophys., 25, 713–729, https://doi.org/10.5194/npg-25-713-2018.
Buehner, M., 2005: Ensemble-derived stationary and flowdependent background-error covariances: Evaluation in aquasi-operational NWP setting. Quart. J. Roy. Meteor. Soc.,
131, 1013–1043, https://doi.org/10.1256/qj.04.15.
——, P. L. Houtekamer, C. Charette, H. Mitchell, and B. He, 2010:
Intercomparison of variational data assimilation and the ensemble Kalman filter for global deterministic NWP. Part II:
One-month experiments with real observations. Mon. Wea.
Rev., 138, 1567–1586, https://doi.org/10.1175/2009MWR3158.1.
Chen, Y., and D. S. Oliver, 2012: Ensemble randomized maximum
likelihood method as an iterative ensemble smoother. Math.
Geosci., 44, 1–26, https://doi.org/10.1007/s11004-011-9376-z.
Cohn, S. E., N. Sivakumaran, and R. Todling, 1994: A fixed-lag
Kalman smoother for retrospective data assimilation. Mon.
Wea. Rev., 122, 2838–2867, https://doi.org/10.1175/15200493(1994)122,2838:AFLKSF.2.0.CO;2.
Courtier, P., J.-N. Thepáut, and A. Hollingsworth, 1994: A strategy
for operational implementation of 4D-Var, using an incremental approach. Quart. J. Roy. Meteor. Soc., 120, 1367–1387,
https://doi.org/10.1002/qj.49712051912.
Emerick, A. A., and A. C. Reynolds, 2012: History matching timelapse seismic data using the ensemble Kalman filter withmultiple data assimilations. Comput. Geosci., 16, 639–659,
https://doi.org/10.1007/s10596-012-9275-5.
——, and ——, 2013: Ensemble smoother with multiple data assimilation. Comput. Geosci., 55, 3–15, https://doi.org/10.1016/
j.cageo.2012.03.011.

Brought to you by University of Maryland, McKeldin Library | Unauthenticated | Downloaded 06/30/21 02:21 AM UTC

2388

MONTHLY WEATHER REVIEW

Errico, R., P. Bauer, and J.-F. Mahfouf, 2007: Issues regarding theassimilation 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.
——, 2018: Analysis of iterative ensemble smoothers for solving
inverse problems. Comput. Geosci., 22, 885–908, https://doi.org/10.1007/s10596-018-9731-y.
——, and P. J. van Leeuwen, 2000: An ensemble Kalman smootherfor nonlinear 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 andtheir implications for data assimilation. Mon. Wea. Rev., 138,
242–255, https://doi.org/10.1175/2009MWR2883.1.
Fairbairn, D., S. R. Pring, A. C. Lorenc, and I. Roulstone, 2014: A
comparison of 4DVar with ensemble data assimilation
methods. Quart. J. Roy. Meteor. Soc., 140, 281–294, https://doi.org/10.1002/qj.2135.
Fletcher, S., 2010: Mixed Gaussian-lognormal four-dimensional
data assimilation. Tellus, 62A, 266–287, https://doi.org/
10.1111/j.1600-0870.2010.00439.x.
——, and A. S. Jones, 2014: Multiplicative and additive incremental variational data assimilation for mixed lognormal
Gaussian errors. Mon. Wea. Rev., 142, 2521–2544, https://doi.org/10.1175/MWR-D-13-00136.1.
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.
Geer, A. J., and P. Bauer, 2011: Observation errors in all-sky data
assimilation. Quart. J. Roy. Meteor. Soc., 137, 2024–2037,
https://doi.org/10.1002/qj.830.
——, and Coauthors, 2017: The growing impact of satellite observations sensitive to humidity, cloud and precipitation.
Quart. J. Roy. Meteor. Soc., 143, 3189–3206, https://doi.org/
10.1002/qj.3172.
——, and Coauthors, 2018: All-sky satellite data assimilation at
operational weather forecasting centres. Quart. J. Roy.
Meteor. Soc., 144, 1191–1217, https://doi.org/10.1002/qj.3202.
——, S. Migliorini, and M. Matricardi, 2019: All-sky assimilation of
infrared radiances sensitive to mid- and upper-troposphericmoisture and cloud. Atmos. Meas. Tech., 12, 4903–4929,
https://doi.org/10.5194/amt-12-4903-2019.
Gelb, A., J. F. Kasper, R. A. Nash, C. F. Price, and A. A.
Sutherland, 1974: Applied Optimal Estimation. The MIT
Press, 374 pp.
Gu, Y., and D. S. Oliver, 2007: An iterative ensemble Kalman filter
for multiphase fluid flow data assimilation. SPE J., 12, 438–
446, https://doi.org/10.2118/108438-PA.
Hamill, T. M., and C. Snyder, 2000: A hybrid ensemble Kalmanfilter-3D variational analysis scheme. Mon. Wea. Rev., 128,
2905–2919, https://doi.org/10.1175/1520-0493(2000)128,2905:
AHEKFV.2.0.CO;2.
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 filter technique. Mon. Wea. Rev.,

VOLUME 149

126, 796–811, https://doi.org/10.1175/1520-0493(1998)126,0796:
DAUAEK.2.0.CO;2.
Jazwinski, A. H., 1970: Stochastic Processes and Filtering Theory.
Academic Press, 376 pp.
Kalman, R. E., 1960: A new approach to linear filtering and prediction problems. Trans. ASME–J. Basic Eng., 82, 35–45,
https://doi.org/10.1115/1.3662552.
Kleist, D. T., and K. Ide, 2015: An OSSE-based evaluation of hybrid variational–ensemble data assimilation for the NCEP
GFS. Part II: 4DEnVar and hybrid variants. Mon. Wea. Rev.,
143, 452–470, https://doi.org/10.1175/MWR-D-13-00350.1.
Liu, C., and Q. Xiao, 2013: Antarctic applications with Advanced
Research WRF using real data. Mon. Wea. Rev., 141, 2721–
2739, https://doi.org/10.1175/MWR-D-12-00130.1.
——, ——, and B. Wang, 2009: An ensemble-based four-dimensional variational data assimilation scheme. Part II: Observing
system simulation experiments with Advanced Research
WRF (ARW). Mon. Wea. Rev., 137, 1687–1704, https://doi.org/10.1175/2008MWR2699.1.
Lorenc, A. C., 2003: The potential of the ensemble Kalman filter
for NWP: A comparison with 4D-Var. Quart. J. Roy. Meteor.
Soc., 129, 3183–3203, https://doi.org/10.1256/qj.02.132.
Lorenz, E. N., 1996: Predictability: A problem partly solved. Proc.
Seminar on Predictability, Reading, United Kingdom, ECMWF,
https://www.ecmwf.int/node/10829.
——, 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/15200469(1998)055,0399:OSFSWO.2.0.CO;2.
Minamide, M., and F. Zhang, 2017: Adaptive observation error
inflation for assimilating all-sky satellite radiance. Mon.
Wea. Rev., 145, 1063–1081, https://doi.org/10.1175/MWRD-16-0257.1.
Morzfeld, M., and D. Hodyss, 2019: Gaussian approximations infilters and smoothers for data assimilation. Tellus, 71A, 1–27,
https://doi.org/10.1080/16000870.2019.1600344.
Neal, R. M., 1996: Sampling from multimodal distributions using tempered transitions. Stat. Comput., 6, 353–366, https://doi.org/10.1007/BF00143556.
Okamoto, K., A. P. McNally, and W. Bell, 2014: Progress towards
the assimilation of all-sky infrared radiances: An evaluation ofcloud effects. Quart. J. Roy. Meteor., 140, 1603–1614, https://doi.org/10.1002/qj.2242.
Poterjoy, J., 2016: A localized particle filter for high-dimensional
nonlinear systems. Mon. Wea. Rev., 144, 59–76, https://doi.org/
10.1175/MWR-D-15-0163.1.
——, and F. Zhang, 2015: Systematic comparison of four-dimensional
data assimilation methods with and without a tangent linear
model using hybrid background error covariance: E4DVar versus
4DEnVar. Mon. Wea. Rev., 143, 1601–1621, https://doi.org/
10.1175/MWR-D-14-00224.1.
——, and J. L. Anderson, 2016: Efficient assimilation of simulated
observations in a high-dimensional geophysical system using a
localized particle filter. Mon. Wea. Rev., 144, 2007–2020,
https://doi.org/10.1175/MWR-D-15-0322.1.
——, F. Zhang, and Y. Weng, 2014: The effects of sampling errors
on the EnKF assimilation of inner-core hurricane observations. Mon. Wea. Rev., 142, 1609–1630, https://doi.org/10.1175/
MWR-D-13-00305.1.
——, R. A. Sobash, and J. L. Anderson, 2017: Convective-scaledata assimilation for the Weather Research and Forecasting
Model using the local particle filter. Mon. Wea. Rev., 145,
1897–1918, https://doi.org/10.1175/MWR-D-16-0298.1.

Brought to you by University of Maryland, McKeldin Library | Unauthenticated | Downloaded 06/30/21 02:21 AM UTC

JULY 2021

KUROSAWA AND POTERJOY——, L. J. Wicker, and M. Buehner, 2019: Progress toward theapplication of a localized particle filter for numerical weather
prediction. Mon. Wea. Rev., 147, 1107–1126, https://doi.org/
10.1175/MWR-D-17-0344.1.
Potter, J., 1964: W matrix augmentation. M.I.T. Instrumentation
Laboratory Memo. SGA, Massachusetts Institute of Technology,
Cambridge, MA, 5–64.
Potthast, R., A. Walter, and A. Rhodin, 2019: A localized adaptive
particle filter 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. Woollen, S. E. Koch, R. Atlas, and R. Hood,
2013: Evaluation of the Earth Systems Research Laboratory’s
global Observing System Simulation Experiment system.
Tellus, 65A, 19011, https://doi.org/10.3402/tellusa.v65i0.19011.
Rommelse, J., 2009: Data assimilation in reservoir management.
Ph.D. thesis, Technical University of Delft.
Sakov, P., D. S. Oliver, and L. Bertino, 2012: An iterative EnKF forstrongly nonlinear systems. Mon. Wea. Rev., 140, 1988–2004,
https://doi.org/10.1175/MWR-D-11-00176.1.
Skjervheim, J.-A., G. Evensen, J. Hove, and J. G. Vabø, 2011: An
ensemble smoother for assisted history matching. SPE
Reservoir Simulation Symp., The Woodlands, TX, SPE, SPE141929-MS, https://doi.org/10.2118/141929-MS.
Snyder, C., T. Bengtsson, P. Bickel, and J. Anderson, 2008: Obstacles
to high-dimensional particle filtering. Mon. Wea. Rev., 136,
4629–4640, https://doi.org/10.1175/2008MWR2529.1.
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.
Thepáut, J.-N., and P. Courtier, 1991: Four-dimensional variational data
assimilation using the adjoint of a multilevel primitive-equation
model. Quart. J. Roy. Meteor. Soc., 117, 1225–1254, https://doi.org/10.1002/qj.49711750206.

2389

van Leeuwen, P. J., and G. Evensen, 1996: Data assimilation andinverse methods in terms of a probabilistic formulation.
Mon. Wea. Rev., 124, 2898–2913, https://doi.org/10.1175/
1520-0493(1996)124,2898:DAAIMI.2.0.CO;2.
Vukicevic, T., T. Greenwald, M. Zupanski, D. Zupanski,
T. V. Haar, and A. S. Jones, 2004: Mesoscale cloud state
estimation from visible and infrared satellite radiances.
Mon. Wea. Rev., 132, 3066–3077, https://doi.org/10.1175/
MWR2837.1.
Whitaker, J., and G. Compo, 2002: An ensemble Kalman smootherfor reanalysis. Proc. Symp. on Observations, Data Assimilation
and Probabilistic Prediction, Orlando, FL, Amer. Meteor.
Soc., 144–147.
——, 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 filter. Mon. Wea. Rev., 132,
1238–1253, https://doi.org/10.1175/1520-0493(2004)132,1238:
IOIEAO.2.0.CO;2.
——, M. Zhang, and J. A. Hansen, 2009: Coupling ensemble
Kalman filter with four-dimensional variational data assimilation. Adv. Atmos. Sci., 26, 1–8, https://doi.org/10.1007/
s00376-009-0001-8.
Zhu, Y., and Coauthors, 2016: All-sky microwave radiance assimilation in NCEP’s GSI analysis system. Mon. Wea. Rev., 144,
4709–4735, https://doi.org/10.1175/MWR-D-15-0445.1.
Zou, X., F. Weng, B. Zhang, L. Lin, Z. Qin, and V. Tallapragada,
2013: Impacts of assimilation of ATMS data in HWRF on
track and intensity forecasts of 2012 four landfall hurricanes.
J. Geophys. Res., 118, 11558–11576, https://doi.org/10.1002/
2013JD020405.

Brought to you by University of Maryland, McKeldin Library | Unauthenticated | Downloaded 06/30/21 02:21 AM UTC