## Abstract

An ensemble of cloud-resolving forecasts from the Weather Research and Forecasting model (WRF) was used to study error covariance for Hurricane Katrina (2005) during a 64-h period in which the storm progressed from a tropical storm to a category-4 hurricane. Spatial error covariance between hypothetical measurements and model state variables was found to be highly anisotropic, variable dependent, and ultimately determined by the underlying storm dynamics, which change dramatically over time. Early in the forecast, when Katrina passed over the southern tip of the Florida Peninsula as a highly asymmetric tropical storm, error covariance structures in the Eulerian coordinates were dominated primarily by position uncertainty, with a secondary dependence on land–air interaction, storm structure, and intensity. The ensemble error dependence on position uncertainty becomes markedly greater with increasing lead time, as diverging storm tracks cause large gradients of wind, temperature, and pressure to be concentrated farther from the mean vortex center.

Ensemble variance for model state variables on storm-relative coordinates becomes increasingly symmetric about the vortex center at greater hurricane intensity. Likewise, spatial and cross-spatial correlations share a similar axisymmetric transition about the origin, while maintaining a large degree of local anisotropy with respect to the location chosen for the correlation point. Our results demonstrate the necessity of using flow-dependent error covariance for initializing a tropical cyclone with dynamically consistent inner-core structure, and provide motivation for future sensitivity experiments pertaining to model resolution and ensemble size.

## 1. Introduction

Improvements have been made over the past three decades in our ability to forecast the track of tropical cyclones (TCs) (see www.nhc.noaa.gov/verification), which is primarily determined by synoptic-scale environmental flow (Wu and Emanuel 1993; Wu and Kurihara 1996; Wang et al. 1998). Despite this progress, our ability to accurately predict intensity change remains quite limited (Elsberry et al. 2007; Houze et al. 2007). The inner-core region of a TC is governed by a multiscale array of dynamic and thermodynamic processes in an air–sea coupled environment. It is apparent that our inability to initialize a TC with dynamically consistent structure and intensity remains a limiting factor for short- to medium-range operational storm prediction. More specifically, small-scale uncertainties associated with moist convection often grow into larger-scale forecast errors in time (Zhang et al. 2003, 2007; Hendricks et al. 2004; Krishnamurti et al. 2005; Zhang 2005; Sippel and Zhang 2008, 2010; Nguyen et al. 2008; Zhang and Sippel 2009). This issue stresses the importance of optimal state estimation, also known as data assimilation, wherein observations, a previous model forecast (prior), and their respective errors are used to estimate a model state prior to integration that contains the least amount of error, denoted the analysis state. Improvements in operational intensity prediction may accompany advancements in tropical cyclone initialization techniques that make better use of in situ and remotely sensed data (Burpee et al. 1996; Leidner et al. 2003; Zhang et al. 2009; Weng et al. 2011).

Data assimilation is usually performed via a *variational* method in which a cost function measuring the distance between a model forecast and observations is minimized, or via an alternative approach in which an *innovation* or correction is added to the prior after it is assigned an optimal weight that minimizes the least squares error (Talagrand 1997). The process of state estimation requires knowledge of both observation and forecast error statistics in terms of variance and covariance. Covariance not only provides an estimate of forecast uncertainty but also quantifies linear multivariate relationships within the model state, thus allowing information to be shared between observed and unobserved variables. In essence, forecast error covariance determines the extent to which a measured quantity can correct state variables in model space. Rather than using stationary, isotropic forecast error covariance, as is the case for standard implementations of variational data assimilation systems, a sequential approach called the ensemble Kalman filter (EnKF) takes into account flow-dependent covariance, typically estimated from an ensemble of short-range forecasts.

Model initialization using an EnKF has gained an increasing amount of attention over the past decade for the analysis and prediction of atmospheric phenomena at regional scales and mesoscales [e.g., Snyder and Zhang 2003; Zhang et al. 2004; Dowell et al. 2004; Tong and Xue 2005; Zhang et al. 2006; Meng and Zhang 2007; Fujita et al. 2007; Meng and Zhang 2008a,b; see Meng and Zhang (2011) for a complete review] and has proven successful for several TC case studies without the use of vortex bogussing (e.g., Chen and Snyder 2007; Torn and Hakim 2009; Zhang et al. 2009; Wu et al. 2010; Weng et al. 2011). In particular, recent experiments in which an EnKF was used to assimilate Doppler radar radial velocity observations into the Weather Research and Forecasting model (WRF) demonstrated a potential to outperform operational prediction systems (Zhang et al. 2009; Weng and Zhang 2011, manuscript submitted to *Mon. Wea. Rev.*).

Bearing in mind the importance of nonstatic, flow-dependent forecast error estimation for dynamically consistent TC initialization, this research is dedicated to the study of state variable covariance for a region that spans the inner and outer core of a developing hurricane. Coherent variance and correlation structures arising from upscale growth of forecast error are examined using an ensemble of forecasts from the WRF for Hurricane Katrina as it intensified from a tropical storm (18 < sustained winds < 33 m s^{−1}) to a category-4 hurricane (59 < sustained winds < 69 m s^{−1}) on the Saffir–Simpson scale. We explore several dynamical features of the cyclones that characterized the ensemble error structures (i.e., vortex asymmetry, tilting, size, intensity, and position uncertainty). In the following section, a brief description of our WRF configuration and method of ensemble generation is provided along with a discussion about our choice of statistics. Results are examined in sections 3 and 4, which include several examples demonstrating how the model dynamics evolve forecast variance and correlations in time. A summary and closing remarks are provided in the last section.

## 2. Forecast model and experiment description

### a. Forecast model and ensemble generation

The fully compressible, nonhydrostatic, mesoscale WRF version 3.1 (Skamarock et al. 2005) was used for this study with a coarse domain of 202 × 181 horizontal grid points at 40.5-km grid spacing and 2 two-way nested domains that automatically follow the storm using the WRF vortex-following algorithm. The innermost domain (D3), where all analysis is performed, has a 253 × 253 horizontal grid with 4.5-km spacing. All domains use 35 vertical levels, most of which are concentrated in the lowest 8 km, with the model top at 10 hPa. The physical parameterization schemes used for this study include the Grell–Devenyi cumulus scheme (Grell and Devenyi 2002) for the two course domains, WRF single-moment six-class microphysics with graupel (Hong et al. 2004), a thermal diffusion scheme for the land surface, the Monin–Obukhov scheme for the surface layer, and the Yonsei State University (YSU) scheme (Noh et al. 2003) to parameterize planetary boundary layer processes.

To generate the first set of model states, ensemble perturbations are produced from the WRF three-dimensional variational data assimilation system (WRFDA) using the cv3 background error covariance option (Barker et al. 2004) and added to the National Centers for Environmental Prediction (NCEP) Global Forecasting System final analysis (FNL) at 0000 UTC 25 August 2005. Lateral boundary conditions are also provided by FNL data and perturbed in the same manner as the initial conditions. The initial 60-member ensemble is integrated for 14.5 h to evolve a flow-dependent forecast error covariance matrix before airborne radar data collected during a flight leg of a National Oceanic and Atmospheric Administration (NOAA) P-3 aircraft are assimilated using an EnKF. Five more legs of radar observations are assimilated at 1530, 1630, 1730, 1900, and 2000 UTC 25 August before an ensemble of 64-h forecasts is produced from EnKF analyses of Hurricane Katrina (Fig. 1). The WRF-EnKF system performed well for Katrina when analyses and the mean forecast are verified against independent flight-level observations and National Hurricane Center (NHC) best-track data; see Weng and Zhang (2011, manuscript submitted to *Mon. Wea. Rev.*) for details regarding the collection and quality control procedures for airborne radar data, as well as the configurations and performance of the EnKF for this event.

The dynamic and thermodynamic structure of forecast members changed dramatically during the simulation as the ensemble mean weakened to a tropical storm between 2000 UTC 25 August and 0000 UTC 26 August before intensifying to a category-4 hurricane by the end of the 64-h forecast. For this study forecast error covariance is evaluated both in a storm-relative reference frame (i.e., relative to a frame of reference following a vortex) and in Eulerian (ground-relative) coordinates. Postprocessing of vortex-following D3 model output was performed by 1) centering WRF data based on minimum sea level pressure (minSLP) and subtracting a storm motion vector from the velocity field of each member for storm-relative calculations, and 2) centering data with respect to a fixed latitude and longitude corresponding to the mean vortex center at the time of the forecast for the Eulerian reference frame calculations. The Eulerian coordinates, used in practice for determining analysis increments in the EnKF algorithm, allow for an investigation of storm position, structure, and intensity effects on forecast error and the resulting EnKF analyses. On the other hand, forecast uncertainty calculated in a frame of reference following the vortex provides further insight into the physical processes forcing the time evolution of the covariance matrix independent of position.

The region of interest for this study lies within the vicinity of the hurricane inner and outer core, where large gradients in pressure, wind velocity, temperature, and moisture contribute to most of the forecast uncertainty. With this in mind, all calculations are performed in a subset of D3 that encompasses the eye, eyewall, and outer rainbands. The new subspace has a grid size of 89 × 89 × 30 and covers a 400 × 400 × 18 km^{3} area, restricting our investigation to a portion of the hurricane where dynamic and thermodynamic relationships are often poorly represented by parameterized covariance in variational data assimilation systems.

### b. Ensemble correlations and variance

Forecast error covariance estimated from ensembles is determined by the underlying error growth dynamics, as discussed by Cohn and Parish (1991), Daley (1992), Evensen (1994), Zhang (2005), and Zhang et al. (2006). When covariance is used for data assimilation, it is informative to separate the statistics into components of correlation and variance. Correlations indicate how information is shared between like and unlike variables in both space and time, while variance quantifies the forecast uncertainty for a given event (variable).

Correlation is a normalized value for covariance and is given by

where *x* and *y* represent two state and/or observed variables at model grid points (*i*, *j*, *k*), and (*i*′, *j*′, *k*′), respectively. Terms denoted with an overbar are ensemble means and *σ* represents sample standard deviations for *x* and *y*. As in Zhang (2005), for the remainder of this paper we will refer to two different definitions of correlation. The first, *cross correlation*, is the element-wise correlation between two different variables at corresponding grid points [*x* ≠ *y*, (*i*, *j*, *k*) = (*i*′, *j*′, *k*′)]. The *spatial* (*auto*-) *correlation* will refer to a correlation between two like variables, calculated with respect to an individual correlation point [*x* = *y*, (*i*, *j*, *k*) ≠ (*i*′, *j*′, *k*′)], and the cross-spatial correlation will be a spatial correlation between two unlike variables [*x* ≠ *y*, (*i*, *j*, *k*) ≠ (*i*′, *j*′, *k*′)].

For data assimilation, the cross-spatial correlation matrix is used to propagate information from an observed scalar quantity to model state variables at neighboring and distant grid points. These statistics combined with uncertainty estimates in terms of observation error and forecast ensemble variance are used to determine the optimal weight associated with two independent approximations of the model state, a forecast state vector denoted by **x*** ^{b}* of length

*N*

**and an observation vector denoted by**

_{x}**y**of length

*N*

**. In the Kalman filter algorithm, the optimal weight is called the Kalman gain**

_{y}**K**and is given by

where is the forecast error covariance matrix of size *N*** _{x}** ×

*N*

**,**

_{x}**H**is an

*N*

**×**

_{y}*N*

**nonlinear observation operator that converts the model state to observation space, and is the observation error covariance matrix (Kalman and Bucy 1961). In the EnKF update equations is parameterized and typically assumed to be a diagonal matrix, and is estimated from an ensemble of**

_{x}*n*forecast states by

A new analysis state vector **x*** ^{a}* is estimated by adding to

**x**

*an increment weighted by*

^{b}**K**:

The analysis increment, defined as the difference between **y** and the forecast state vector in observation space **Hx*** ^{b}*, determines the direction in which a state variable is corrected (i.e., a decrease or increase in value), whereas

**K**determines the amplitude. In the traditional implementation of the EnKF,

**x**

*and*

^{a}**x**

*in (4) are taken as analysis and forecast members.*

^{b}The impact of forecast spatial covariance on (2) and (4) can be more easily understood by considering the case where a single observation is assimilated. For a scalar observation *y*′ with variance *R*′, from (2) reduces to an *N*** _{x}** × 1 column matrix , given by cov(

**x**

*,*

^{b}**x**

^{b}**H**

^{T}), and reduces to a scalar

*d*. The EnKF update equation now becomes

where corrections to the prior forecast now depend on and the difference between *y*′ and the model forecast of *y*′, given by **Hx*** ^{b}* (Snyder and Zhang 2003).

Several types of EnKF algorithms exist that deviate in one or more ways from the above equations, such as the ensemble square root filter version of Whitaker and Hamill (2002) used for this study, which requires additional calculations for updating the ensemble with unique estimates of **K** for each member. Nevertheless, the fundamental idea provided by (2)–(4) and the single observation case in (5) is adequate for exemplifying the significance of background (forecast) error covariance for the success of data assimilation.

## 3. Error covariance during the first 16 h

We first analyze the first 16 h of forecast uncertainty (along with the EnKF analysis ensemble at 2000 UTC 25 August), using both the Eulerian and storm-relative reference frames. In this section, we compare error covariance from the analysis with that of the 4-, 10-, and 16-h forecasts (0000–1200 UTC 26 August). The analysis error at 2000 UTC 25 August represents the posterior uncertainty provided by the EnKF after assimilating several hours of inner-core airborne Doppler radar velocities for Hurricane Katrina [see details in Weng and Zhang (2011), manuscript submitted to *Mon. Wea. Rev.*]. Most members followed a southwestward track after initialization and decreased in intensity upon making landfall over the southern tip of the Florida peninsula. In general, members were in large agreement with Katrina’s intensity change throughout the first 10 h of model integration.

Using a vortex-following coordinate system (with the location of minSLP for each member relocated to the origin), Figs. 2a–d show that the standard deviations *σ _{υ}* of meridional wind velocity

*υ*decreased by about 4 m s

^{−1}within the radius of maximum wind (RMW), corresponding to the decrease in intensity during this period (Fig. 1c). Despite the reduction in forecast uncertainty within the RMW,

*σ*contours of 4 m s

_{υ}^{−1}or greater increased spatially at the lowest model level (~35 m in our experiment) with greater forecast lead time. Forecast uncertainty in the low-level wind field was most likely caused by moist convection—in particular, nonlinear interactions associated with surface sensible and latent heat fluxes, which were aided by the spread of storm positions. While the mean maximum 10-m wind speed

*V*

_{max}decreased by 6 m s

^{−1}from the analysis to 4-h forecast, the ensemble mean minSLP remained nearly constant because of a lag in the mass field response to the falling wind speeds. As a result, standard deviations

*σ*of storm-relative pressure

_{P}*P*decreased slightly at the lowest model level (Figs. 2e–h) and remained less than 6 hPa throughout the first 10 h of integration. Likewise, large temperature gradients in the inner-core region, extending outward from the eye, concentrated the largest ensemble temperature

*T*uncertainty near the mean vortex center, as indicated by

*T*standard deviations

*σ*in Figs. 2i–l. Maximum

_{T}*σ*values were found near 1.9 km in altitude and decreased substantially by 0600 UTC 26 August as most members tracked near land, decreasing the upward flux of moist air from the ocean surface, latent heat release in the eyewall, and horizontal temperature gradients at middle levels. With position error removed in the storm-relative calculations, it is apparent that error associated with intensity and structural uncertainty of Katrina’s inner-core circulation decreased throughout the 10-h forecast. Nevertheless, as a majority of members moved farther from the coast of southwest Florida by 1200 UTC 26 August, ensemble

_{T}*σ*and

_{P}*σ*increased substantially.

_{T}Ensemble analysis and forecast *σ _{υ}* are again shown in Fig. 3, but with all calculations performed in the original (Eulerian) model coordinates. In this case, analysis error evolves from a centralized, compact structure to a broad, asymmetric distribution after a short forecast period. Analysis

*υ*uncertainty is maximized inside the RMW, with

*σ*increasing inward from 4 to 20 m s

_{υ}^{−1}in the marine boundary layer (MBL) because the initial distribution of storm positions concentrated large wind gradients near the mean vortex center. Given the small analysis position spread, forecast uncertainty initially resembles that of the storm-relative reference frame (Fig. 2a) but rapidly grows within the first several hours of model integration as members diverge from the mean position. Unlike the storm-relative calculations, the spatial extent of 4 m s

^{−1}or larger forecast

*σ*in Eulerian coordinates grew substantially with increasing lead time, and an uncertainty maximum of approximately 24 m s

_{υ}^{−1}formed near the mean storm center by 16 h. A notable characteristic of low-level

*σ*is an approximately 2 m s

_{υ}^{−1}decrease between sea surface and land grid points. A majority of ensemble members passed over the southern tip of Florida as tropical storms and winds over land were damped by surface roughness, resulting in a coastal variance discontinuity for an approximately 150-m-deep layer (three WRF

*η*levels in our experiments). This is consistent with Zhang et al. (1999), who observed a sharp ocean-to-land decrease in surface winds observed from their numerical simulation of Hurricane Andrew (1992).

Similar to *σ _{υ}*, uncertainty with regard to Katrina’s simulated thermodynamic fields (

*P*and

*T*) in the Eulerian coordinates rapidly broadens with increasing lead time and track spread (Fig. 4). A clear depiction of modes within the storm track can be observed in three-dimensional distributions of

*σ*and

_{P}*σ*(not shown), since the largest anomalies for these variables typically occur near storm centers instead of along two or more regions in the RMW (as is usually the case for

_{T}*υ*). In general, the error structures associated with

*P*and

*T*are highly sensitive to non-Gaussian distributions of storm position, which becomes increasingly noticeable with growing track spread. To isolate the effects of position error, a “track error only” (TEO) ensemble was created by calculating a storm-relative mean from the 60-member WRF ensemble forecast and then azimuthally averaging the mean and placing 60 identical wavenumber-0 vortices in the original locations of members valid at the forecast time. The TEO ensemble presents a hypothetical case in which all members are axisymmetric and share the same structure; therefore, only position uncertainty exists. Figure 5 compares

*P*analysis and forecast statistics throughout the first 16 h for the Eulerian and TEO ensemble experiments. From this comparison it is evident that most of the asymmetric structure in

*σ*can be attributed to track spread, as the storm-relative error acts only to enhance the magnitude of

_{P}*σ*in the Eulerian reference frame. Despite the analysis position uncertainty being mostly Gaussian, the limited sample size for this experiment allows the eastward located outliers to create a secondary

_{P}*σ*maximum that appears in both Eulerian and TEO calculations (Figs. 5a,e). By 0000 UTC 26 August, three

_{P}*σ*maxima have formed (Figs. 5b,f), corresponding to three independent storm position modes with errors of 4 hPa or greater near the surface. The three maxima persist after 10 h of integration (Figs. 5c,g), but two predominant modes of storm track appear, one on each side of the Florida Keys corresponding to

_{P}*σ*> 10 hPa by 1200 UTC 26 August (Figs. 5d,h).

_{P}Besides the use of variances in statistical data assimilation systems for determining the amount of confidence associated with a given state estimation, an observation’s impact on correcting a prior forecast also depends critically on the cross and/or spatial correlations between observed and unobserved variables [(5)]. Multivariate spatial correlations that illustrate the effects of assimilating hypothetical observations to correct 10-h surface pressure and zonal winds are chosen as examples. When vortex position in terms of central latitude and longitude is chosen as the observed value to be assimilated (Figs. 6a,e), large correlation dipoles form, with values that alternate sign across the domain. Assuming a bivariate normal distribution for ensemble positions, the eigenvectors and eigenvalues accounting for 69% (*σ* = 65.8 km) and 31% (*σ* = 44.2 km) of the position standard deviations are plotted in Figs. 6a and 6e, with the leading component oriented along the direction of largest variance (primary axis of the distribution). The dipole resulting from corr(central longitude, *P*) lies in nearly the same direction as the leading eigenvector (the primary axis), since eastward located members tend to be positioned farther northward and westward members tend to be more southward. The subsequent correlations would act to increment the location of minSLP for Katrina members forward or backward along the primary axis of the position distribution for this case, rather than along lines of constant longitude. Similarly, corr(central latitude, *P*) would increment members in the direction of the secondary axis. The result is consistent with the analysis increments described in Chen and Snyder (2007), where vortex position was assimilated using an EnKF with a barotropic model, except pressure is used here as the variable to be updated rather than vorticity.

In Fig. 6b, a zonal wind velocity *u* located at a grid point approximately 75 km west of the mean vortex center at 0600 UTC 26 August (denoted observation point a and measurement *u*_{a}) is cross-correlated spatially with the *P* field, yielding a broad positive signal in the lower-left portion of the domain directly south of point a. If the WRF-EnKF analysis/forecast cycle were to continue from the 10-h forecast, by (5) measurement *u*_{a} would produce positive (negative) analysis increments south of point a for an observation greater than (less than) the mean forecasted *u*_{a}, thus filling (deepening) Katrina members with centers near the given location. At the same observation point, corr(*u*_{a}, *u*) (Fig. 6f) are physically consistent with the hypothetical *u*_{a}*-*to*-P* corrections, with positive values north of the largest *P* increments. Consequently, a negative pressure increment from *u*_{a} would coincide with corrections to the velocity field that increase the cyclonic rotation around the largest *u*_{a}*-*to*-P* signal. Similarly, for a measure of *u* ~ 30 km southwest of the mean storm center (point b), corr(*u*_{b}, *P*) in Fig. 6c indicate a strong negative signal north of the mean vortex center. Despite its proximity to location a, measurement *u*_{b} provides a completely new set of information regarding the forecast *P* and *u* fields that would result in EnKF corrections to the model state in regions of the domain that are dissimilar to those from measurement *u*_{a}, thus demonstrating the tremendous amount of variability that exists in the forecast covariance matrix. Once again, the increments produced from *u*_{b} will act to increase (decrease) the cyclonic rotation near the pressure falls (rises) provided by the spatial correlations shown in Figs. 6c and 6g. In a third example, a wind observation taken at location c, a distance greater than 100 km from the storm center (Figs. 6d,h), would yield corrections to the *P* field on both sides of the mean vortex center, much like the case in which storm position is assimilated. Recall that error covariance is determined by both multivariate spatial correlations and forecast uncertainty, and consequently the EnKF increments resulting from the provided measurements would also depend on *σ _{P}* and

*σ*in these examples.

_{υ}As demonstrated in Fig. 5, the spatial distribution of ensemble variance depends greatly on storm position uncertainty. To extend our investigation into the effects of track error on ensemble forecast error covariance, the TEO ensemble was again used to calculate the same spatial and cross-spatial correlations shown in Fig. 6. The resulting signals (cf. Fig. 7) are remarkably similar to the pure Eulerian case, implying that the assimilation of storm position, along with point measurements *u*_{a}, *u*_{b}, and *u*_{c}, acts mostly to correct vortex location. Central latitude and longitude observations force the members along the eigenvectors shown in Figs. 7a and 7e and similar *P* increments result from wind observations, only the axis in which the correlation dipole forms is dependent on the location of the wind measurement with respect to the mean storm center. The spatial correlations calculated for point measurements of *u* will act to move the prior vortex center in a manner consistent with the central latitude and longitude correlations, provided that a sufficient amount of observations are available. As more observations are assimilated, the influence of position uncertainty on the ensemble background (forecast) covariance will decrease, as the storms are incremented closer to the center realized by the observations. The amount of complexity associated with the covariance matrix, even for an ensemble of identical wavenumber-0 vortices, is demonstrated by the distinct spatial correlation structures that exist for the hypothetical measurements.

Despite the similarities between the Eulerian and TEO ensemble correlations, it is evident that forecast uncertainty with regard to intensity and structure still has a significant effect on error covariance, even at lead times as large as 10 h. For example, storms that are positioned farther north by 0600 UTC 26 August encounter more land–air interaction and experience less of a relationship between *u* and distant *P* values. As a result, the TEO ensemble overestimates the magnitude of negative correlations in the northeastern quadrant of the domain in Fig. 7b, where *u*_{a} only has a small bearing on *P* anomalies over this region for the pure Eulerian case. Similar arguments can be applied for explaining cases in which the TEO ensemble underestimates signals in the correlation matrix as well, such as in Fig. 7c, where track error alone cannot account for the large linearity that exists between decreasing *u*_{b} and pressure rises on the southwestern coast of Florida (cf. Fig. 6c).

## 4. Storm-relative forecast error covariance beyond 16 h

It is apparent from the previous section and past studies analyzing EnKF analysis increments (i.e., Chen and Snyder 2007; Torn and Hakim 2009) that model forecast error is primarily determined by position spread, especially for large lead times. At the same time, error covariance resulting from uncertainty with regard to hurricane intensity and dynamical structure is still quite significant and must be accounted for during state estimation. Despite the practical importance of calculating ensemble error statistics in an Eulerian reference frame for vortex-scale data assimilation, errors emerging from sources other than position uncertainty are considerably more difficult to visualize in the Eulerian coordinates. For that reason, a storm-relative reference frame will be used for the remainder of this paper to study forecast error covariance past 16 h of model integration. Forecast times representative of major changes in Katrina’s wind field and thermodynamic error structures were chosen for the storm-relative analysis, corresponding to mean storm intensity forecasts of category-1 (1200 UTC 26 August), category-2 (0000 UTC 27 August), category-3 (1200 UTC 27 August), and category-4 hurricanes (1200 UTC 28 August).

While viewing spatial correlations calculated with respect to a storm-relative reference frame, the reader should be cognizant of the fact that the correlation point chosen for the calculation is in reality not fixed in the ground-relative sense. In other words, the ground-relative position of a hypothetical observation to be assimilated will differ for each member. Despite this limitation for ground-relative data assimilation purposes, calculations of storm-relative forecast error covariance have a large functional importance for scenarios in which observations are assimilated in short time intervals where vortex position uncertainty is minimal, such as the 1-h analysis/forecast cycles used for assimilating radar observations in the beginning of our simulation. They also provide a cleaner comparison between flow-dependent and parameterized forecast covariance typically used in variational data assimilation methods that make use of a vortex bogussing algorithm for correcting vortex position before any observations are assimilated.

### a. Wind and temperature

As indicated by the dashed lines in Figs. 1b,c, intensity error for the Katrina forecast ensemble grew steadily throughout the 64-h forecast (with the exception of the initial weakening due to interaction with land throughout the first 10 h), as standard deviations of *V*_{max} and minSLP increased from 5.7 to 8.5 m s^{−1} and 7.6 to 24.6 hPa, respectively. Following the initial 0–10-h decrease in maximum storm-relative wind variance (Figs. 2a–c), *σ _{υ}* continued to grow while organizing into two maxima on opposite ends of the eyewall near corresponding maxima in wind speed. In association with low-level radial inflow that becomes more axisymmetric about the origin with hurricane intensity, the axis separating positive and negative values of mean

*υ*(black contours in Fig. 8) contains a varying degree of tilt throughout the forecast. This can be conceptualized by separating the

*υ*field into tangential and radial components and realizing that the radial contribution is positive in the southern half of the domain and negative in the northern half, while the tangential contribution has positive and negative maximums directly east and west of the origin, respectively. In Figs. 8a–d, vectors representing the radial contribution to the mean

*υ*field (from this point forward denoted

*υ*) are plotted to illustrate how inflow at this level affects the mean

_{r}*υ*field and distribution of

*σ*. For most of the members, a significant amount of southerly inflow persists throughout the simulation. The inflow contributes to the asymmetric shape of the mean

_{υ}*υ*field in the southern half of the domain and explains the region of

*σ*> 6 m s

_{υ}^{−1}that extends southwestward from the origin. It is apparent from these forecast times that regions of large forecast

*υ*uncertainty (

*σ*> 6 m s

_{υ}^{−1}) existing outside the RMW tend to coincide with regions where

*υ*is at a maximum. As the

_{r}*υ*vectors increase in magnitude for the northern half of the domain by 1200 UTC 28, the radial inflow within 100 km of the origin becomes largely axisymmetric, causing the orientation of the mean

_{r}*υ*field north of the origin to match that of the southern portion of the domain. The ensemble-estimated forecast uncertainty shows a large sensitivity of

*σ*to asymmetries in radial inflow, thus again highlighting the value of flow-dependent covariance modeling in vortex-scale data assimilation.

_{υ}The vertical structure of *σ _{υ}* also exhibited considerable changes throughout the course of the simulation. When most Katrina members were near tropical-storm and category-1 hurricane intensity (1200 UTC 26 August–0000 UTC 27 August),

*σ*> 6 m s

_{υ}^{−1}occupied a volume bounded by the eyewall (~45 km from vortex center) extending as high as 17 km in altitude (Figs. 8a,b,e,f). With minimal intensity spread at these times, wind error may have resulted from high wavenumber asymmetries caused by convection and/or vortex tilting, as weaker members are prone to synoptic-scale environmental forcing and vertical wind shear. By 1200 UTC 27 August (Figs. 8c,g), most of the members obtained category-3 intensity and became increasingly axisymmetric. The increasing dependence on wavenumber-0 vortex structure led to weaker

*σ*in the eye but larger

_{υ}*σ*near the RMW as the ensemble intensity forecasts diverged and uncertainty associated with the magnitude and location of maximum wind speeds became more prevalent. The transition is clearer in Figs. 8i–l, where correlations between minSLP and the

_{υ}*υ*field are shown. At 1200 UTC 26 August and 0000 UTC 27 August, the portions of the wind field varying with minSLP are loosely organized and largely asymmetric. Significant correlations of magnitude greater that 0.6 are absent from the inner 27 km of the vortex center at mid- to upper levels because of a lack of linearity between intensity and the wind field at these heights (perhaps due to vortex tilting). It is also possible that land–air interaction early in the forecast acted to diminish low-level correlations at 1200 UTC 26 August. Once the ensemble mean surpassed category-2 intensity, the resulting standard deviations and correlations organized into symmetric structures about the origin, which had only minor fluctuations in magnitude and spatial extent between 1200 UTC 27 August and 1200 UTC 28 August.

Unlike the variance observed for wind components, the distribution of forecast temperature uncertainty remained mostly axisymmetric throughout the entire forecast, only increasing in magnitude along large thermal gradients in the eyewall. Figure 9 shows maximum *σ _{T}* increasing from 2.6 to 3.2 K in regions of the eye between 1200 UTC 26 August and 0000 UTC 27 August. Early in the simulations, most of the error was confined to lower and middle levels (

*z*< 5 km), but adiabatic warming induced by the thermally indirect secondary circulation raised the

*σ*maximum with greater lead times. The large column of

_{T}*σ*> 4.2 K, extending from 4 to 9 km at 1200 UTC 28 August, was aided by the intensity spread at later times, since the presence of both weak and strong members produced varying degrees of warming in the eye. A secondary uncertainty maximum at 14 km emerged between 1200 UTC 27 August and 1200 UTC 28 August as a result of enhanced upper-level warming for strong members at later times when the simulated Katrina members obtained a mean intensity greater than category 3.

_{T}As demonstrated in Figs. 8 and 9, ensemble forecast uncertainty induced by changes in Katrina’s intensity and structure are greatest in the vicinity of the eyewall and eye region near large horizontal and vertical gradients of wind, pressure, and temperature. Since the purpose of data assimilation is to make corrections to the model state vector, given that an observation is correlated with the system’s state variables, hypothetical observations will again be chosen to analyze the value certain measurements hold for improving the EnKF analysis. Our objective in this section is to highlight important aspects of the system-scale vortex uncertainty that may produce a large signal in the Eulerian reference frame but are often difficult to distinguish in the covariance matrix given the presence of position uncertainty.

Figure 10 shows the significance of assimilating temperature observations in the eye of Katrina at altitudes of 6.3 (location d) and 3.4 km (location e). The *T* spatial correlations at location d (*T*_{d}) show midlevel warming being positively correlated with *υ* in eastern portions of the domain and negatively correlated with *υ* in western regions with respect to the origin. The same pattern emerges in Figs. 8a–h in the form of two horizontally tilted *υ* error maxima. It is also significant to mention that the three-dimensional distributions of corr(minSLP, *υ*) (only the vertical profile is shown in Figs. 8i–l) produces a qualitatively similar signal due to the hydrostatic relationship between upper-level warming and low-level pressure falls in the eye. For the low-level measure of *T* at location e (*T*_{e}) (Figs. 10i–l), the signal is similar to that provided from *T*_{d} for the earlier time steps, albeit weaker, and is nearly nonexistent by 1200 UTC 28 August. This suggests that warming of the eye at 3.4 km is indicative of stronger members during early stages of the storm’s intensification but has less importance at later times where the mean reaches category-4 intensity. From Fig. 11, where a 14-km temperature at location f is correlated with the *υ* field, it is clear that by 1200 UTC 28 August the wind fields of the simulated hurricanes are more associated with upper-level warming than at lower and middle levels. The correlations resulting from the three temperature values show that the height at which this coupling occurs is largely dependent on the intensity of the simulated hurricanes. They also reveal that a significant amount of information can be spread from middle- and upper-level measures of temperature for correcting the magnitude and structure of the wind field at nearly all layers of the model. Flow-dependent covariance calculated from these hypothetical measurements would lead to corrections in Katrina’s wind field that are dynamically consistent with the model’s uncertainty with regard to changes in storm structure and intensity. The resulting EnKF increments would be highly anisotropic and take into account warming trends in the hurricane eye, as well as the primary and secondary circulation of the wind field, given the tilted covariance structures associated with inflow.

Figure 12 demonstrates another hypothetical scenario in which low-level tangential velocity *V _{θ}*, a wind velocity derivable from satellite or airborne Doppler radar, is used for correcting the wind field (denoted measurement

*V*

_{θ}_{g}). Spatial correlations with respect to

*V*at a point approximately 50 km from the vortex center (denoted location g) extend from the lower to middle levels at 16 h, with significant values confined to the eastern half of the domain inside a radius of 100 km. From 1200 UTC 26 August to 1200 UTC 28 August, correlations of magnitude greater than 0.4 extended outward from the center, becoming increasingly axisymmetric with larger lead times. The signal grew considerably in magnitude, with correlations of 0.8 or greater persisting in the RMW for levels as high as 12 km by 1200 UTC 27 August. By 1200 UTC 28 August, the RMW of most members contracted to a distance inside the hypothetical storm-relative location marked by point g, thus causing a slight decrease in correlations at this time; that is, the location of point g with respect to the RMW is changing throughout the forecast, which also influences the magnitude and structure of the correlations in Fig. 12.

_{θ}Using the same surface wind value, cross-spatial correlations with respect to *V _{θ}*

_{g}and the

*T*field are shown in Fig. 13, yielding a complex relationship between storm intensity and the inner- and outer-core thermodynamic structure of the simulated Katrina. In general, the correlations are most significant near the surface and in the storm’s eye region, with strong anticorrelations appearing in regions outside the eyewall. Figures 13a–d show a clear relationship between surface wind speeds and low-level temperatures in an extensive portion of the domain throughout the simulations. As surface air from outside the eyewall circulates inward from regions in the northern portions of the domain, the more intense Katrina members benefit from the inflow of warmer, more buoyant air from these regions, resulting in positive correlations. In the middle and upper portions of the domain, a column of positive correlations inside the eyewall results from members with larger sustained surface wind speeds coinciding with a strong secondary circulation and increased subsidence and warming in the eye. The vertical distribution of these correlations is similar in shape to the

*T*standard deviations shown in Fig. 9, producing the most significant amount of covariance for this example.

Anticorrelations of notable magnitude exist for a large, isolated region directly outside the eyewall and below the freezing level throughout the simulations. A consistent yet slightly stronger signal appears for correlations between minSLP and the *T* field (not shown), demonstrating that cooler temperatures in this portion of the domain coincide with stronger storms. While most of the deep convection and large concentrations of rain, snow, ice, and graupel occur upwind of this region, the advection and melting/evaporation of precipitation into the northwestern quadrant of the outer core suggest one possible explanation of the observed correlations. It is also feasible that melting or evaporation of hydrometeors occurred before reaching the area of large correlations, allowing cooler air resulting from the phase change to be advected into this region. As the storm obtains a more axisymmetric profile by 1200 UTC 28 August, the large negative signal vanishes, possibly due to the larger updrafts surrounding northwestern portions of the eyewall by this time (the vertical velocity field will be discussed in the following section). Considering the low ensemble variance collocated with these correlations (*σ _{T}* < 1 K) (cf. Fig. 9), the amount of cooling found outside the eyewall in this region for strong members is relatively insignificant; therefore, the correlations would yield only marginal updates during the data assimilation procedure.

Since significant correlations occupy regions of large forecast variance in this example, temperature error can be effectively corrected using measurement *V _{θ}*

_{g}in a storm-relative reference frame. The resulting EnKF increments would be most significant inside the hurricane eye as opposed to a region surrounding point g inside the eyewall, which would not be the case for a traditional variational data assimilation system using isotropic forecast error covariance. Examples such as this, along with the largely asymmetric wind/temperature variance and correlation fields at early forecast times (Figs. 8–13), further validate the importance of flow-dependent forecast error covariance for adequate initialization of developing hurricanes.

### b. Moisture and vertical motion

Diabatic energy sources in terms of sensible and latent heat fluxes in the MBL are important for the steady-state maintenance and intensification of TCs (Rotunno and Emanuel 1987). Furthermore, intensity and vortex core asymmetries can be highly sensitive to the distribution of tropospheric moisture (Nolan et al. 2007; Nguyen et al. 2008; Sippel and Zhang 2008, 2010; Fang and Zhang 2010, 2011). Under the above realizations, an examination of ensemble forecast uncertainty of water vapor mixing ratio *q _{υ}* and vertical velocity

*w*is well justified.

The distribution of storm-relative, low-level *q _{υ}* (Fig. 14) shows a large amount of moist air in the southeastern quadrant at earlier times that becomes increasingly axisymmetric throughout the forecast. Radial gradients of

*q*change sign and increase substantially in the eye, with a maximum directly inside the RMW aided mostly by 1) the vertical transfer of moisture by updrafts in the eyewall, 2) evaporation of hydrometeors at dry midlevels inside the RMW, and 3) an upper-level minimum in the eye region, where subsidence of dry upper-level air occurs. By 1200 UTC 26 August (Figs. 14a,e) most of the

_{υ}*q*forecast uncertainty

_{υ}*σ*

_{qυ}is found in the MBL near areas of deep convection, such as the eyewall and spiral rainbands. As ensemble intensity forecasts diverge with increasing lead time, a larger uncertainty emerges with regard to the upward flux of moist air in the eyewall and dry subsidence in the eye. As a result,

*σ*

_{qυ}inside the RMW continued to grow, surpassing 2.4 g kg

^{−1}by the end of the forecast, while standard deviations in rainbands outside the eyewall remained nearly constant.

Evidence of a moisture uncertainty dependence on storm intensity is shown in Fig. 15, where low-level *q _{υ}* in the eye and upper-level

*q*in the eyewall is negatively correlated with minSLP. The strongest negative signal in the correlation matrix appears in a 4-km-deep layer above the 12-km level, where stronger hurricanes transport moister air upward and out from the origin at all time steps, resulting in positive

_{υ}*q*anomalies above the eyewall that spiral outward with height. Likewise, a relationship between surface pressure falls and middle- to upper-level drying in the eye forms at later time steps as the Katrina members intensify and form well-developed warm cores. The height at which drying in the eye is correlated with minSLP moves upward with increased storm intensity, which is consistent with Figs. 10 and 11, which show a height dependence on storm intensity and warming. It is not prudent to interpret these correlations as evidence of a warm core height dependence on intensity. They simply suggest that as hurricanes in the WRF intensify, stronger members tend to coincide with anomalously warmer and dryer air at higher altitudes.

_{υ}Near the surface, distributions of corr(minSLP, *q _{υ}*) show a low-level moisture dependence on deepening pressure and intensity. Significant correlations outside the eyewall are limited to discrete regions of the domain throughout the simulations, with the exception of the last time step at 1200 UTC 28 August when the Katrina ensemble reached its peak intensity. As members strengthen to category-4 intensity, a larger signal emerges between low-level

*q*and minSLP because of the increased evaporation and a more organized transport of moist air toward the center. This causes an increase in moist convection and heating that is consistent with the observed pressure falls at the forecast time. An increase in radial inflow in the northern section of the domain may have also contributed to the larger correlations, as moist air in this region is more efficiently transported inward by the end of the simulation (cf. Fig. 8).

_{υ}As mentioned earlier, the inner core is characterized by intense slantwise updrafts and convection in the eyewall. At larger radii, bands of organized convection associated with regions of positively buoyant air and upward motion exist alongside clear areas of downdrafts, most likely associated with the return circulation from outflow. From 1200 UTC 26 August to 1200 UTC 28 August, the original asymmetric, spiraling updrafts (and downdrafts) near the vortex center (Figs. 16a,e) gave way to a ring of strong *w* about a calm eye (Figs. 16d,h). As Katrina members intensified and became increasingly axisymmetric, forecast uncertainty with regard to *w* in the eyewall decreased. At earlier stages of the simulation, the location of maximum updrafts with respect to the vortex center differed substantially among members (recall that wavenumber-0 vortices were centered based on minSLP for storm-relative calculations), resulting in standard deviations of *w* (*σ _{w}*) reaching a maximum by 1200 UTC 26 (Figs. 16a,e). Outside the eyewall, the magnitude of

*σ*became saturated (

_{w}*σ*< 2.0 m s

_{w}^{−1}) near discrete convective cells associated with the outer rainbands (Figs. 16b,c,d,f,g,h).

Despite the importance of *w* in the initialization of convective-scale features, the correlations between vertical motion and other variables are small in our experiment, even at earlier time steps (although this may not have been the case if the ensemble were initialized at a time in which Katrina was a major hurricane). The only variable that correlated considerably with *w* for lead times greater than 16 h (in the storm-relative case) was *w* itself (Fig. 17). In addition, only grid points inside strong updrafts within the eyewall were able to generate strong enough signals to warrant substantial corrections to the *w* field. Since vertical motion is tied closely to transformations of hydrometeor species and their effects on buoyancy, correlations in which one or both variables are rain, snow, ice, or graupel mixing ratio are also small in this forecast ensemble (not shown). The loss of correlation between *w* or hydrometeor species with other state variables is likely because *w* and hydrometeor species have a more significant amount of small-scale variability and are more directly associated with moist convection, which is largely nonlinear and inherently less predictable (Zhang et al. 2006). One possible means of extracting significant vortex-scale *w* cross correlations from the ensemble, such as those associated with the secondary circulation, is to remove small-scale features from *w* or hydrometeor species via spectral analysis or other filtering tools before estimating the correlations. This additional step of postprocessing model output, which will be investigated by the authors in a follow-up study, may allow for more dynamically consistent corrections to be made to the low-wavenumber (more balanced) component of the vortex structure of the *w* field.

## 5. Summary and conclusions

In an effort to improve our current understanding of ensemble forecast covariance pertaining to TCs, correlations and standard deviations of model state variables were analyzed from a 60-member WRF ensemble of Hurricane Katrina during a 64-h period in which forecasts progressed from a tropical storm to a category-4 hurricane. The examination was restricted to a 400 × 400 × 18 km^{3} model subspace that captured the vortex-scale variance and correlations. Previous research (e.g., Houze et al. 2007; Nguyen et al. 2008; Fang and Zhang 2011) suggests that the dynamical elements of this region may have a profound influence on intensity fluctuation, intensification, and decay. Our study was motivated by the small- to large-scale error propagation that occurs when inner-core features of an initialized TC vortex are structurally inconsistent or ill accounted for during model analysis time steps.

The observed error statistics were shown to be dependent on ensemble track spread, land–air interaction, storm structure, and intensity. Position spread was quickly realized as the largest factor for determining the observed flow-dependent variance and correlation distributions, even at relatively short forecast lead times (e.g., 4 h). Members were also centered based on the location of minSLP for storm-relative calculations to separate the effects of track spread and analyze secondary sources of forecast uncertainty. At longer lead times (corresponding to larger intensity), the resulting storm-relative covariance was less affected by high wavenumber asymmetries and had a larger dependence on intensity uncertainty.

Our results demonstrate how a better understanding of TC error growth dynamics for complex regions, such as the inner core, can be achieved through the analysis of model forecast ensembles. The resulting spatial covariance matrices, determined by the underlying model dynamics, are increasingly axisymmetric about the storm center for storm-relative calculations as Katrina progresses into the mature hurricane stages but remain largely anisotropic with respect to the correlation points. It appears that sufficient correlations exist between standard measured variables and the wind, temperature, and pressure fields in locations of high forecast uncertainty; thus, adequate flow-dependent updates can be made to the given state variables during the data assimilation procedure. Correlations involving moist processes and buoyancy-driven vertical velocity are typically lower but still significant at the vortex scale. The flux of heat and moisture from the ocean surface to the lower troposphere, accompanied by discrete convective cells outside the eyewall and water transformations in the middle and upper troposphere, tends to complicate possible links between these variables and other predictors. Results from both ground- and storm-relative calculations suggest that ensemble-based data assimilation methods can adequately regulate inner-core structure and intensity, provided that track spread is marginal or a suitable amount of data has been previously assimilated to decrease forecast position uncertainty.

While the scope of this study exclusively covers the change in ensemble covariance for a developing hurricane, the results allow for speculation into the covariance for nonintensifying storms. Calculations made from the Katrina ensemble result mostly from balanced vortex dynamics, with additional flow-dependent, asymmetric structure caused by forcing from smaller-scale features. Since most of the forecast uncertainty and multivariate correlations can be attributed to the primary and secondary circulation and to the adiabatic structure of the developing hurricanes, large similarities are expected for steady-state or weakening tropical cyclones, where similar dynamics apply. Nevertheless, one should always be cautious with results obtained from a single-case study. In future research, the authors will extend their investigation to more storms including steady-state or weakening cases. The Katrina ensemble provides valuable information with regard to the targeting of observations through spatial/cross-spatial correlations, demonstrates the necessity of short assimilation/analysis cycles for adequate intensity and structure corrections, and verifies the advantage of implementing high spatial resolution data for vortex-scale data assimilation. This study presents a descriptive look at how ensemble-based data assimilation methods extract statistical information from a sample of forecasts to make dynamically and thermodynamically consistent corrections to a hurricane simulation. It highlights advantages and potential shortfalls of the ensemble Kalman filtering method and its application to tropical cyclone data assimilation.

Future work will involve more direct comparisons between analysis and short-term forecast error covariance and sensitivity with respect to model resolution and ensemble size. The increase in resolution may (or may not) lead to improvements for small-scale vertical velocity and moisture correlations, as highly nonlinear error growth mechanisms, such as moist convection, may limit such correlations to shorter ensemble lead times. Likewise, a larger ensemble size may be required in sync with the increased model resolution to take full advantage of the higher degrees of freedom associated with each forecast. Given the complexity of the forecast error covariance matrix, as demonstrated in many of our examples, future studies using idealized hurricane models will be necessary to fully understand how model dynamics and physical parameterization schemes force the observed ensemble error propagation.

## Acknowledgments

The authors thank Yonghui Weng for performing the Katrina ensemble simulations for this study and Daniel Stern for his helpful comments and proofreading. This work was supported in part by the NOAA Hurricane Forecast Improvement Project (HFIP), Office of Naval Research Grant N000140910526, and the National Science Foundation Grant ATM-0840651. The computing was performed at the Texas Advanced Computing Center.