Two four-dimensional hybrid data assimilation systems based on the Weather Research and Forecasting (WRF) Model are applied over the life cycle of Hurricane Karl (2010). One method uses a mix of ensemble- and climate-based background error covariance in a four-dimensional variational data assimilation (4DVar) system that uses an adjoint model to assimilate observations over a time window (denoted E4DVar). The second method approximates the function of linearized models in 4DVar with perturbations generated from an ensemble forecast using the full nonlinear model (denoted 4DEnVar). Ensemble perturbations in 4DEnVar provide a four-dimensional covariance, which is combined with a static climate-based covariance for performing the data assimilation. In cycling data assimilation experiments, analyses produced by both methods provide more accurate intensity forecasts than E3DVar, owing mostly to a better representation of moisture near the developing storm. Despite providing a computationally efficient alternative to E4DVar, predictions made from 4DEnVar analyses are less accurate than E4DVar for the tested case study. Numerical experiments using identical background error statistics in both schemes reveal differences in the mesoscale structure of the developing storm, which are suspected to be responsible for this result. In particular, 4DEnVar analyses contain a less intense inner-core circulation and lower column relative humidity than E4DVar at analysis times closest to Karl’s genesis, which lead to a persistent slow bias in intensifying the storm. These results suggest errors introduced in the linearization of the model for E4DVar may be less severe than errors introduced by the localization of time covariances. This study provides the first comparison of hybrid E4DVar and 4DEnVar for a tropical cyclone application.
Data assimilation approaches that combine strategies from variational and ensemble Kalman filtering (EnKF) methods are applied regularly for operational weather forecasting and research. Applications include global numerical weather prediction (e.g., Buehner 2005; Buehner et al. 2010a,b; Bishop and Hodyss 2011; Clayton et al. 2013; Kuhl et al. 2013; Wang et al. 2013; Wang and Lei 2014; Lorenc et al. 2015) as well as regional analysis and forecasting (e.g., Liu et al. 2009; Zhang and Zhang 2012; F. Zhang et al. 2013; Liu and Xiao 2013; Poterjoy and Zhang 2014a, hereafter PZ14a; Gustafsson and Bojarova 2014). One common result from these studies is that variational data assimilation methods benefit from using either an ensemble-based estimate of prior error covariance in place of the climate-based static covariance, or from a hybrid covariance that combines ensemble and static background errors. This finding has been verified by systematic comparisons of hybrid and nonhybrid data assimilation systems as well (e.g., F. Zhang et al. 2013).
Among the ensemble-variational data assimilation strategies listed above are four-dimensional (4D) methods, which seek the maximum likelihood model solution, given a prior estimate of the model state or background, a set of observations over a time window, and the uncertainty of both quantities. The first operational use of 4D data assimilation for weather prediction came in the form of incremental four-dimensional variational data assimilation (4DVar; Courtier et al. 1994; Rabier et al. 2000), which applies tangent linear and adjoint versions of the forecast model to propagate information in time to assimilate observations displaced over a window. This approach uses a model forecast from the previous analysis as the background state, and approximates the background uncertainty using a climate-based error covariance (Parrish and Derber 1992). The tangent linear model propagates increments, or differences between the control variable and a background state, along a nonlinear trajectory to future times, while the adjoint model propagates sensitivity gradients backward from observation times to the beginning of the window. The linear model operators evolve the background error covariance implicitly in time, thus providing a source of flow-dependent information for the state estimation (Lorenc 2003a).
When ensemble information is used in 4D data assimilation systems, a choice emerges regarding how to model the temporal covariances required for performing the state estimation. The EnKF, first proposed by Evensen (1994), presents a means of propagating error covariances forward in time without linearized model operators. EnKFs rely on an ensemble of model forecasts to approximate the prior mean and error covariance, which are updated each time observations are available. Algorithms that expand the EnKF to include a time dimension include the ensemble Kalman smoother (Evensen and van Leeuwen 2000) and 4DEnKF (Hunt et al. 2004). Likewise, ensemble 3DVar data assimilation systems have been conceived by Lorenc (2003b), Zupanski (2005), and Buehner (2005). Extending these methods to 4D, Liu et al. (2008) introduced the 4D-ensemble-variational (4DEnVar) technique, which uses an ensemble of model trajectories to find the maximum likelihood solution. 4DEnVar contains a number of qualities that make it a practical alternative to 4DVar as a data assimilation method for atmospheric models. Minimizing the 4DEnVar cost function does not require the coding of tangent linear and adjoint models, which is appealing for operational centers and research institutions that do not currently maintain a 4DVar system. An equally important factor is the scalability of the algorithm on massively parallel computing platforms; unlike 4DVar, 4DEnVar does not require serial integrations of a tangent linear model and its adjoint during each iteration of the cost function minimization.
Despite the advantages 4DEnVar presents for implementation, methods that use ensemble covariance in 4DVar are found to be more accurate for global numerical weather prediction models (Buehner et al. 2010b; Lorenc et al. 2015). One challenge for 4DEnVar comes from the localization of 4D ensemble covariances; this problem is avoided in the 4DVar framework because a full-rank localized covariance is evolved implicitly in time by linearized models.1 The lack of linear model operators in 4DEnVar also complicates the process of evolving errors from the climate-based static error covariance when a hybrid background error covariance is used. How to model this covariance component in an efficient manner remains an open research question at the time of writing this manuscript. Using a hybrid covariance in 4DVar, Poterjoy and Zhang (2015) find that the ability of the tangent linear and adjoint model to evolve temporal information from a static covariance component is most advantageous for treating model errors not accounted for by the ensemble. The static covariance improves the analysis in a manner similar to additive covariance inflation in ensemble filters: It provides additional orthogonal directions for errors to grow in the assimilation window, but with the added advantage of remaining full rank. Nevertheless, 4DEnVar may contain benefits for certain weather regimes that have not been realized in past global modeling applications. For example, the linearized model in 4DVar often contains simplified physics to avoid strong nonlinearities in parameterization schemes (Sun and Crook 1997; Honda et al. 2005; Huang et al. 2009), which is avoided in 4DEnVar. Evidence of 4DEnVar producing more accurate analyses than ensemble implementations of 4DVar for regional-scale modeling have been demonstrated by Gustafsson and Bojarova (2014). At the same time, the 4D ensemble introduces nonlinearities in the temporal evolution of covariances, possibly leading to a poorly conditioned minimization problem (J. Sun 2014, personal communication).
In a recent study by Poterjoy and Zhang (2014b), the authors apply a coupled EnKF–4DVar (E4DVar) data assimilation system to assimilate routinely collected observations and field measurements taken during the development and intensification of Hurricane Karl (2010). Measurements assimilated during this study were collected during the Pre-Depression Investigation of Cloud-systems in the Tropics (PREDICT) field campaign. Poterjoy and Zhang (2014b) provide a detailed introduction into this application by discussing the spatial and temporal frequency of observations, along with their impact on synoptic and mesoscale processes leading to Karl’s genesis and intensification in the Weather Research and Forecasting (WRF) Model. Poterjoy and Zhang (2014b) also highlight the important role of mesoscale processes in the evolving weather event, which motivated PZ14a to test a newly developed E4DVar system for generating analyses of Karl’s entire life cycle.
The E4DVar method used in PZ14a incorporates a hybrid ensemble/climate-based background error covariance in a 4DVar system developed for the WRF Model (Zhang et al. 2014). Each E4DVar data assimilation cycle includes a hybrid update of the ensemble mean state using 4DVar, and an EnKF update of the ensemble perturbations. The authors find E4DVar to produce a more accurate representation of wind, temperature, and moisture near the pre-Karl disturbance than similarly configured 4DVar and EnKF data assimilation systems, which leads to better predictions for the timing and location of genesis and rapid intensification. In the current study, the authors apply a coupled EnKF–4DEnVar system for the same application to examine the effectiveness of a 4D-ensemble data assimilation system that does not require linearized models. Results are compared with a coupled EnKF–3DVar (E3DVar) data assimilation system and the E4DVar method presented in PZ14a. The purpose of this study is to examine E3DVar, E4DVar, and 4DEnVar analyses for an application where theoretical differences highlighted in Poterjoy and Zhang (2015) may have a large impact on the modeling of 4D covariances by E4DVar and 4DEnVar. For this reason, the authors focus on a single well-observed weather event. The goal of this evaluation is to supplement future and existing studies that provide more rigorous statistical comparisons of adjoint- and ensemble-based 4D data assimilation methods.
The organization of the manuscript is as follows. We describe the model, data assimilation systems, and cycling data assimilation setup in section 2. In section 3 we present a set of single-observation experiments using ensemble- and climate-based background error covariance for observations displaced at multiple times in the assimilation window. Sections 4 and 5 contain results from cycling and noncycling data assimilation experiments, respectively, where noncycling refers to the use of the same background at each analysis time. We provide our conclusions and a discussion in section 6.
2. Experiment setup
a. WRF Model setup and case study
This study uses version 3.4.1 of the Advanced Research WRF Model (Skamarock et al. 2008). The model domain for ensemble forecasts has 251 × 226 horizontal grid points with a spacing of 13.5 km (see Fig. 1 of PZ14a) and 35 vertical levels with a 5-mb (1 mb = 1 hPa) upper boundary. Deterministic forecasts initialized from each analysis use a 4.5-km spacing two-way nested domain of size 253 × 253 that follows the storm with preset moves. We use the following physical parameterization schemes for unresolved processes in the model: Monin–Obukhov similarity (Monin and Obukhov 1954) for the surface layer, the Yonsei University planetary boundary layer scheme (Noh et al. 2003), five-layer thermal diffusion for surface layer physics, the Rapid Radiative Transfer Model (RRTM; Mlawer et al. 1997) and Dudhia (Dudhia 1989) radiation schemes, and the WRF single-moment 6-class microphysics scheme (Hong et al. 2004). Both domains represent cumulus convection explicitly, which was found in offline sensitivity experiments to provide more accurate temperature and moisture fields than simulations using parameterized convection for this model setup and case study.
Each data assimilation method tested during this study uses ensemble forecast data to estimate background errors. The initial ensemble of model states is formed at 1800 UTC 7 September by adding random perturbations generated from a month-long estimated climatological background error covariance to the National Centers for Environmental Prediction (NCEP) Global Data Assimilation System (GDAS) analysis at this time. Following a 12-h ensemble forecast, the first data assimilation cycle occurs at 0600 UTC 8 September, and observations are assimilated every 6 h. The data assimilation systems are cycled over a period that includes the days leading up to Karl’s development at 1200 UTC 14 September, as well as its initial intensification on 15 September. Deterministic forecasts from each cycle are run to 0000 UTC 18 September, which spans the entire life cycle of Karl. Observations assimilated during these experiments include regularly collected surface and upper-air data from the National Oceanic and Atmospheric Administration (NOAA) Meteorological Assimilation Data Ingest System (MADIS) and dropsonde measurements taken during the PREDICT field campaign (Montgomery et al. 2012). Observations from the National Aeronautics and Space Administration (NASA) Genesis and Rapid Intensification Processes (GRIP) campaign (Braun et al. 2013) are combined with MADIS and PREDICT soundings to verify deterministic forecasts from each method. See Poterjoy and Zhang (2014b) and Torn and Cook (2013) for details regarding the impact of the PREDICT field observations on analyses and forecasts for this event.
b. Hybrid 3DVar, 4DVar, and 4DEnVar
For the variational component of the data assimilation, we use version 3.4.1 of the WRF data assimilation (WRFDA) package (Huang et al. 2009; Barker et al. 2012). For a given observation time t, the variational analysis comes from minimizing a cost function that measures the misfit of a control variable to the background state and observations . The minimization is carried out with respect to increments from (Courtier et al. 1994), and the cost function is expressed as the sum of background () and observation () terms. The incremental 3DVar cost function depends on quantities valid at the data assimilation time, represented by time 0:
where and are the background and observation error covariance matrices, respectively. The quantity is the innovation vector, which is given by the difference between and the background projected into observation space by the observation operator . Likewise, increments are projected into observation space using the tangent linear version of , denoted by . Equation (1) is extended to 4DVar for observations displaced temporally between and :
In (2), the control variable is the vector of increments at the beginning of an assimilation window , which relates to increments at future times through the tangent linear forecast model .2 The time displacement of observations is also reflected in the 4DVar innovation vector, which is given by
where the background state at is propagated forward to the time of observations using the nonlinear forecast model . See X. Zhang et al. (2013) and Zhang et al. (2014) for detailed descriptions of the WRF 4DVar data assimilation system and the WRF tangent linear model, respectively.
To include a hybrid covariance in the variational cost function, is separated into two terms:
The first and second terms of (4) are the climatological and ensemble contribution to the analysis increment, respectively (Wang et al. 2008b). Each in the second term is a Schur product between the vector and the nth ensemble perturbation, where . Following Lorenc (2003b), the ensemble information is introduced in the variational cost function using an additional term for the ensemble background information ():
The α control variables in the second term on the right-hand side of (5) store the weighting coefficients () for each ensemble perturbation, which are constrained by ; matrix is block diagonal, and holds the correlations used to localize the ensemble covariance. The climate and ensemble background cost terms are weighted by the scalar coefficients and to determine the contribution of the two background error covariances during the minimization. The quantities and are chosen to satisfy , which guarantees the weights assigned to the climate- and ensemble-based covariance sum to unity3 (Wang et al. 2008b). Equation (5) provides the cost function used in E4DVar experiments for introducing ensemble information in 4DVar, which reduces to E3DVar for τ = 0.
To avoid the use of a tangent linear model, (6) can be approximated as
where is neglected from the climate contribution (first term on the right-hand side), and the propagation of ensemble perturbations by (second term on the right-hand side) is replaced with perturbations from a nonlinear ensemble forecast. The assumptions made in the first and second terms in (7) are equivalent to those made by 3DVar first guess at appropriate time (FGAT; Fisher and Andersson 2001) and 4DEnVar (Liu et al. 2008, 2009), respectively. The hybrid 4DEnVar system applied in this study minimizes (5) using the approximations made in (7) to incorporate static and ensemble covariance in the cost function. Unlike Liu et al. (2008, 2009), ensemble perturbations are stored in model space and projected into observation space using the tangent linear observation operator.
c. Coupled data assimilation
The three data assimilation methods tested in this study require the execution of both EnKF and variational data assimilation systems—which can be performed simultaneously. From an ensemble forecast, the members are separated into mean and perturbations, which are introduced into the variational system as the background state and perturbation vectors , respectively. To estimate the analysis state at the center of the window for E4DVar and 4DEnVar, we propagate the sum of analysis increments and background state forward from the beginning of the window.4 Assuming all observations displaced through the assimilation window are valid at the analysis time (center of the window), the EnKF estimates the posterior perturbations, which are added to the maximum likelihood solution found during the variational analysis. The posterior ensemble is then used to initialize the ensemble forecast for the next data assimilation cycle. This approach results in a two-way coupling of ensemble and variational data assimilation methods (Zhang et al. 2009) for E3DVar, E4DVar, and 4DEnVar.
d. Data assimilation configuration
For the 4D data assimilation experiments, observations are binned into 1-h intervals and assimilated over a 6-h period that spans ±3 h from the analysis time. Likewise, the E3DVar experiment assimilates all observations over the same period, but assumes all data are valid at the same time. We also applied E3DVar with 3-h assimilation cycles (and a 3-h window), but results from this experiment are omitted from the manuscript because no significant improvements were found over the 6-h configuration. All three methods use the same climatological background error covariance, which is estimated over the previous month using 24- and 12-h forecast differences (Parrish and Derber 1992). We use the default covariance amplitudes and length scales for the climatological errors, since tuning these parameters yielded no significant change in performance. The ensemble contains 60 members, and its covariance is localized in the variational and EnKF components using a 900-km horizontal radius of influence and 34 vertical levels for vertical localization. Here, the variational methods perform a model-space Gaussian localization using a recursive filter (Wang et al. 2008a), while the EnKF applies a Gaspari and Cohn (1999) localization function directly in the Kalman gain matrix when assimilating each observation.
Ensemble spread between data assimilation cycles is maintained by inflating the posterior ensemble perturbations after each analysis using the “covariance relaxation” technique described in Zhang et al. (2004). This method uses 80% of the prior perturbations and 20% of the posterior perturbations to update the posterior ensemble members. The data assimilation sensitivity to localization and inflation parameters was explored in offline experiments before deciding on the current configuration. Our choice of weighting coefficients for the climate and ensemble background cost terms come from results of sensitivity tests performed in previous regional modeling studies (Zhang and Zhang 2012; F. Zhang et al. 2013). The specified weights translate into a hybrid configuration that uses 20% of the static covariance and 80% of the ensemble covariance. Extensive tests of hybrid E4DVar and 4DEnVar in Poterjoy and Zhang (2015) using a low-dimensional dynamical system show the methods having a similar sensitivity to covariance localization and inflation, window length, ensemble size, and hybrid weighting coefficients. As a result, the optimal set of parameters for the two methods were found to be near identical for almost all model and observing system scenarios tested. While the parameter sensitivity of both methods may not be as comparable for complex applications, like tropical cyclogenesis in WRF, this result provides a basis for using the same parameter values for E4DVar and 4DEnVar, which is done in the current study. Maintaining the same parameter values between experiments also allows for a simpler interpretation of results. As applied in this study, the only difference between these two methods is the modeling of covariances via the tangent linear and adjoint operators or an ensemble of model trajectories.
To reduce the computational cost of running the tangent linear and adjoint models for the 4DVar component of E4DVar, we perform the analyses using a reduced 40.5-km resolution to minimize the cost function (Zhang et al. 2014). We also generate the background trajectory for calculating innovations and linearizing the model and observation operators only once for each cycle (i.e., we do not run multiple outer loops). Because the 4DEnVar and E3DVar analysis steps require less computing time, we make use of the full-resolution 13.5-km grid for the minimization in these experiments.
3. Single-observation experiments
E4DVar and 4DEnVar differ primarily in how each method models background error covariance in the assimilation window. In this section, we explore these differences by assimilating synthetic observations and comparing the resulting analysis increments. We apply E4DVar and 4DEnVar using ensemble and static covariance to assimilate single observations displaced at three different times in the assimilation window. The comparison is performed for the 1200 UTC 13 September cycling time using the background (ensemble forecast) from the cycling E4DVar experiment. This cycle occurs 24 h prior to genesis, during a time in which the low-level mesoscale vortex undergoes a large increase in circulation (see discussion in section 4a). We create 850-mb temperature observations at 0, 3, and 6 h from the analysis time by taking temperature values from a 6-h forecast from the beginning of the assimilation window and adding a 1-K error. We then assimilate the observations separately using E4DVar and 4DEnVar (Fig. 1) to show potential temperature (θ) increments at the beginning of the assimilation window associated with observations at times t = 0, t =3, and t = 6 h. These experiments are performed using 100% and 0% contributions of ensemble covariance to demonstrate how each method estimates time covariances from the ensemble and static covariance components.
E4DVar and 4DEnVar provide identical analysis increments for observations at the beginning of the assimilation window (t = 0 h, Figs. 1a,f), but clear differences emerge for observations placed at future times. The θ increments change sign over the window, becoming negative for observations at t = 3 h and t = 6 h. We suspect that uncertainty in vortex position contributes the most to this sign change: negative increments produced from assimilating observations later in the assimilation window contain increasingly larger scales, which serve to correct track error inferred from future observations. Nevertheless, intensity errors are also large at this stage of the storm, which further complicates the interpretation of these experiments. The increments in the top row of Fig. 1, are produced using ensemble covariance only. In this case, the differences between increments for t > 0 h result from the following: 1) the unique localization strategies applied by the methods—E4DVar localizes the ensemble covariance at the beginning of the window, then propagates information forward using linearized models, whereas 4DEnVar performs the localization after estimating time covariances; and 2) simplifications introduced in the WRF tangent linear and adjoint models used for E4DVar. For the application chosen to perform these comparisons, the effects of using a different localization strategy are accelerated by the rapid development of the storm over the assimilation window. The sensitivity of E4DVar and 4DEnVar to forecast error growth was explored previously in Poterjoy and Zhang (2015) using the 40-variable Lorenz (1996) system with varying degrees of model forcing. Results from this study show that faster forecast error growth in the assimilation window leads to a faster degradation of 4DEnVar analyses, because of the localization strategy used to remove spurious temporal correlations. This finding explains some of the differences seen in the 100% ensemble covariance cases plotted in the top row of Fig. 1. Likewise, the elevated forecast error growth rate for this application may also accelerate the effects of errors introduced by imperfect linearized models in 4DVar, which may degrade the accuracy of E4DVar increments in this example.
In the other limiting case, where 0% of the ensemble covariance is used, increments plotted in the bottom row of Fig. 1 reflect the contribution of the static background error covariance. The E4DVar increments in this row (Figs. 1f,g,i) are identical to traditional implementations of 4DVar using a climate-based background covariance, while the 4DEnVar increments (Figs. 1f,h,j) are identical to 3DVar-FGAT. As the observation is displaced further in time, the structure of the E4DVar increments becomes closer to the 100% ensemble covariance case (i.e., a large area of negative θ increments forms around the developing storm). The result follows from the tangent linear model and its adjoint having more time to evolve flow-dependent information over the assimilation window. Comparisons using hybrid covariance in the two systems are omitted from Fig. 1, because the resulting increments are a linear combination of the 100% and 0% ensemble cases described in this section [cf. (4)].
4. Results from cycling experiments
a. Analysis storm structure
In this section, we compare E3DVar, E4DVar, and 4DEnVar analyses over a weeklong cycling period between 0600 UTC 8 September and 0000 UTC 15 September, which spans the development and intensification of Karl. For this comparison, we focus on the kinematic and thermodynamic characteristics of analyses in the vicinity of the storm. As described in PZ14a, we estimate the pregenesis storm center by averaging the 950- and 700-mb circulation centers, which are designated as the positions that maximize the azimuthal mean winds within a 3° radius at each level. If a tropical storm is present in the simulation (10-m winds exceeding 18 m s−1), we perform a Barnes’s analysis (Barnes 1964) using the 10-m, 850-mb, and 700-mb vorticity fields to locate the storm center. Figures 2 and 3 summarize the evolution of the pregenesis disturbance in each set of analyses by comparing the time series of area-averaged quantities from 1200 UTC 9 September to 0000 UTC 15 September, following a short spinup period between 0600 UTC 8 September and 0600 UTC 9 September.
Figures 2a,b compare 950- and 500-mb relative vorticity (ζ) values averaged within 3° of the storm center. In all three experiments, the evolution of the low- and midlevel circulation at the meso-α scale is qualitatively consistent with the synopsis provided by Stewart (2010), likely because of the large spatial and temporal extent of in situ observations assimilated. The low-level circulation increases on 11 September, followed by a decrease between 11 and 13 September before strengthening until Karl intensified into a tropical storm at 1800 UTC 14 September (indicated by the vertical dashed line in Fig. 2). Meanwhile, the strength of the midlevel vortex increases steadily throughout the entire cycling period (Fig. 2b). The development of the 950-mb vortex between 13 and 15 September in the analyses coincides with the alignment of the low- and midlevel circulation centers; this alignment is indicated by the time series of the 950–500-mb tilt magnitude plotted in Fig. 2c. The vortex alignment in the analyses leads to a reduction in the local 950–500-mb vertical wind shear (Fig. 2d), which we estimate from winds within 3° of the storm center. The plotted shear magnitude takes into account the environmental shear, as well as the effects of asymmetries induced by the tilted vortex. Mechanisms that can contribute to the alignment of a vortex in vertical wind shear are discussed in Schecter et al. (2002), Nolan and McGauley (2012), Rappin and Nolan (2012), and Zhang and Tao (2013). Poterjoy and Zhang (2014b) found vortex alignment to be an important factor in determining when tropical cyclogenesis occurs in simulations initialized from EnKF and 4DVar analyses of this event, and noted the importance of dropsondes near the pregenesis disturbance in capturing the low-level development of the vortex with time.
With the exception of a faster vortex alignment in the E4DVar analyses, the three data assimilation methods yield qualitatively similar kinematic structure for meso-α-scale features of the disturbance preceding genesis. One notable distinction between analyses produced from the 4D data assimilation methods and E3DVar is the amount of near-saturated air surrounding the developing storm. Figure 2e compares the time series of spatially averaged 950–500-mb column relative humidity (CRH), which is defined as the ratio of vertically integrated water vapor to vertically integrated saturation water vapor. Over the cycling period examined here, the E3DVar analyses remain 2%–8% (RH) drier than E4DVar and 4DEnVar, which inhibits the generation of convectively induced vorticity anomalies in the pregenesis disturbance, and decreases the likelihood of genesis (Hendricks et al. 2004). The relatively low values of CRH near the developing storm in E3DVar analyses are comparable to EnKF and 4DVar analyses discussed in Poterjoy and Zhang (2014b) for Karl (see their Fig. 3). In addition, verifications of deterministic forecasts in section 4b show that E4DVar and 4DEnVar improve the representation of thermodynamic variables over E3DVar. These results suggest that the hybrid 4D data assimilation methods capture moist processes related to Karl’s genesis more accurately than the ensemble filters and non-ensemble methods tested for this case study (PZ14a; Poterjoy and Zhang 2014b).
Figure 3 compares the same quantities plotted in Fig. 2, except values are averaged within 1° of the storm center. Here, distinctions emerge with regard to the inner core of the storm in each set of analyses. The 13–15 September assimilation cycles reveal a gradual intensification of the midlevel vortex in the E4DVar and 4DEnVar data assimilation experiments (Fig. 3b); the E3DVar analyses contain a similar increase in circulation at these times, except the vortex remains much weaker. Following the increase in circulation in the middle troposphere, the low-level vortex intensifies rapidly between 14 September and the end of the cycling period on 15 September. The intensification of the inner-core midlevel vortex in E4DVar and 4DEnVar analyses follows a 24-h period of elevated CRH values. These magnitudes agree with the 80% average relative humidity threshold found in WRF simulations by Nolan (2007) to be a precursor for genesis. Furthermore, the cycle at which the initial intensification occurs in each experiment is consistent with the rate at which the atmosphere approaches saturation in the volume used to perform these calculations. The E4DVar experiment produces elevated values of CRH earlier in the cycling period and develops a mesoscale vortex faster than E3DVar and 4DEnVar in the analyses.
b. Deterministic forecasts
Here we compare deterministic forecasts initialized with the analyses described in section 4a. The prediction accuracy in each experiment is quantified by comparing forecasts with routine soundings and 234 total dropsonde measurements collected during PREDICT and GRIP. We again focus on meteorological fields in the vicinity of the storm by calculating forecast root-mean-square differences (RMSDs) to observations within 800 km of the National Hurricane Center (NHC) best-track center.5 This verification is performed for zonal and meridional wind (u and υ), temperature (T), and water vapor mixing ratio (). We average the RMSDs over the volume that spans the verification region between 950 and 200 mb, and temporally over 28 cycles between 0600 UTC 8 September and 0000 UTC 15 September. Results are plotted for forecast times between 0 and 72 h in Fig. 4, with vertical bars that indicate the 99% confidence interval, assuming the samples of RMSDs are uncorrelated in time and distributed normally. Figure 5 provides a more detailed overview of forecast errors for 24- and 72-h lead times by showing forecast bias and RMSDs averaged over 9 vertical columns between 1000 and 200 mb (observations within 50 mb of each level are used). This figure also indicates the number of verifying observations for each level on the right-hand side of each plot.
The 0-h verification in Fig. 4 measures the fit of observations to the analysis. Because the verifying observations include data assimilated by all three schemes (i.e., measurements from MADIS and PREDICT), a component of the mean RMSDs at the initial time comes from a model state that is dependent on the verifying observations. For this reason, we omit an interpretation of the analysis accuracy using the 0-h RMSDs. The time series of errors beginning from 0 h, however, provides some evidence regarding imbalance in the initial conditions for the two 4D methods, which require a 3-h integration from the beginning of the assimilation window. For example, the T RMSDs decrease immediately after 0 h in the 4DEnVar experiments. Likewise, 4DEnVar produces the largest 0-h RMSDs for u, υ, and q, which we verified to be caused in part by elevated gravity wave activity in the 3-h forecasts to the analysis time (not shown). Previous studies using 4DEnVar in global models have encountered similar issues related to balance because of the use of a discontinuous set of data to characterize the 4D prior statistics (e.g., Kleist and Ide 2015). Lorenc et al. (2015) and Buehner et al. (2015) both found benefits using four-dimensional incremental analysis update (4DIAU) to control noise during initialization. The E4DVar system used in this study uses a weak digital filter constraint 4DVar that penalizes high-frequency oscillations in the cost function, but no such mechanism is available in the WRF-4DEnVar system used in this study. PZ14a used the digital filter constraint in their comparison of 4DVar and E4DVar because of minor benefits exhibited by 4DVar with the additional cost function term. This benefit occurs mostly because the increments produced during inner-loop iterations depend largely on a climatological covariance estimate, which represents background errors rather poorly for tropical cyclone applications. These increments result in an elevated level of imbalance, which can be dealt with using the digital filter. Because this problem is partially alleviated by the improved covariance estimates in E4DVar, the role of this filter becomes less important. Offline tests using E4DVar without the digital filter constraint show no significant benefit of this additional cost function term, so we use the same E4DVar experiment in PZ14a and the current paper to maintain consistency between the two studies. While repeating the E4DVar experiment without the digital filter would provide a more consistent comparison between the methods, we have no reason to suspect the digital filter changes the major results of this study.
We use short-term forecasts after 0 h to verify the effectiveness of the data assimilation methods in producing accurate predictions of state variables in the region surrounding the disturbance. From this verification, the 4DEnVar and E3DVar experiments have comparable predictive skill for u, υ, and T out to 72 h. Forecasts from E4DVar analyses contain similar υ RMSDs over the same forecast range, but have significantly lower u and T RMSDs than E3DVar and 4DEnVar for lead times greater than 48 h. The lower u and T RMSDs in E4DVar forecasts are found mostly in the upper model levels (Fig. 5), and coincide with lower biases than E3DVar and 4DEnVar. On the other hand, E4DVar produces the largest bias for υ in low levels, despite having RMSDs that are comparable to E3DVar and 4DEnVar for most model levels. Both 4D data assimilation methods produce a more accurate representation of moisture than E3DVar near the disturbance, which is reflected in the significantly lower RMSDs for starting from 12 h. Figures 5g,h indicate that most of the forecast error in these experiments is due to a negative bias, which is smaller for the E4DVar and 4DEnVar cases at most levels. This result suggests that the elevated values of CRH found near the disturbance in E4DVar and 4DEnVar analyses, compared to E3DVar in Figs. 2e and 3e, provide a better depiction of Karl’s pregenesis relative humidity field.
We also verify the predicted track and intensity of the developing cyclone in deterministic forecasts. Here, the performance of the data assimilation techniques depends on how well each method captures features in analyses that produce accurate predictions of the storm track before and after genesis, as well as the correct timing of genesis and rapid intensification in the simulations. Figures 6 and 7 show deterministic track and intensity (maximum 10-m wind speed) forecasts, initialized every 6 h between 1800 UTC 12 September and 0000 UTC 15 September. Simulations are run through the entire life cycle of Karl, thus providing insight into the long-range forecast performance of the three data assimilation methods. Black lines plotted in Figs. 6 and 7 represent the observed track and intensity, which are taken from the NHC best-track dataset; pregenesis storm positions come from the wave-tracking product described in Wang et al. (2012). Table 1 also provides the mean RMSEs, calculated from all deterministic forecasts plotted in these figures at 0-, 6-, 12-, 24-, 48-, and 96-h forecast lead times.
For forecasts initialized 3–4 days before genesis (Fig. 6), the E4DVar forecasts propagate the pregenesis wave westward, and produce accurate genesis predictions as early as 1200 UTC 12 September, 54 h before Karl became a tropical storm. Forecasts from 4DEnVar and E3DVar analyses fail to capture the correct translation speed and direction of the wave until 1800 UTC 11 September, which coincides with the second PREDICT flight mission. Likewise, none of the 4DEnVar and E3DVar forecasts prior to 13 September capture Karl’s genesis before making landfall on the Yucatan Peninsula. Figure 7 compares the track and intensity forecasts from the three experiments within two days of Karl becoming a tropical storm. During this period, all three experiments produce a persistent right-of-track forecast bias, which also occurred in operational hurricane model forecasts for this event (Stewart 2010). The variability in Karl’s track between data assimilation cycles after 1800 UTC 11 September also remains relatively low, which is consistent with operational Global Forecast System forecasts—used in this study for lateral boundary conditions. Nevertheless, the E4DVar forecast errors are consistently smaller than forecast errors from 4DEnVar and E3DVar analyses. In addition to providing the most accurate track forecasts, the E4DVar analyses produce superior forecasts for the timing of Karl’s development, which leads to better predictions of the storm’s rapid intensification over the Bay of Campeche. The faster development of Karl in E4DVar analyses—compared to 4DEnVar and E3DVar—is demonstrated by the 1° averaged vorticity comparisons in Figs. 3a,b. In the 4DEnVar and E3DVar simulations, Karl’s development is delayed until the weather system crosses the Yucatan Peninsula, which leads to a persistent low bias in the intensity forecasts (Fig. 7).
5. Results from noncycling experiments
For the experiments analyzed in section 4, the mean and covariance from a given data assimilation cycle is propagated forward to the next cycle and used to define the new prior error distribution. The cycling causes errors to accumulate between analyses because of deficiencies in how each method estimates posterior quantities. One negative outcome of the cumulative errors is that each method may respond differently to random noise introduced by sampling errors in the prior ensemble statistics and observation errors, thus affecting the outcome of the results presented in section 4. The impact of cumulative errors during cycling can be especially large for tropical cyclogenesis cases, where the underlying dynamics are very sensitive to small-scale initial condition errors in kinematic and thermodynamic variables such as vorticity and humidity (e.g., Sippel and Zhang 2008, 2010; Zhang and Tao 2013; Munsell et al. 2013). To remove differences due to cumulative errors introduced by the data assimilation system and observation noise, we perform a new set of analyses and forecasts using the same prior ensemble in the three data assimilation systems. For these experiments, each method uses forecast members from the cycling E4DVar experiment discussed in section 4. We generate new analyses and deterministic forecasts for 10 cycles between 1800 UTC 12 September and 0000 UTC 15 September, which covers 48 h before Karl’s genesis, as well as its initial intensification on 15 September. Results from these experiments depend exclusively on each method’s unique data assimilation strategy, which allows for an investigation into possible systematic differences in the three methods.
We also performed experiments using a hybrid configuration of 4DEnVar that does not include the climate-based background error covariance in the hybrid cost function (denoted 4DEnVar-ep). Instead, we introduce by replacing the ensemble forecast perturbations at the beginning of the window with hybrid perturbations:
where each is drawn randomly from This approach allows 4DEnVar to evolve perturbations from the static contribution of the background error covariance in the assimilation window to approximate the time evolution of Poterjoy and Zhang (2015) demonstrate advantages of this application of 4DEnVar for cases when substantial model error is present in the forecast system. Representing with hybrid perturbations has a similar effect as additive inflation because it allows forecast errors to grow in new orthogonal directions during the assimilation window. One disadvantage, however, is that the method introduces additional sampling error during the data assimilation because of the replacement of a full rank with an ensemble representation.
Deterministic track and intensity forecasts initialized from each of the four hybrid analyses are plotted in Fig. 8, along with forecasts from the background mean state used to perform each analysis (denoted control throughout this section). When no data are assimilated to improve the background state, the forecasts generate persistently accurate track forecasts, but contain a slow bias in the intensification of Karl (Figs. 8a,b). Forecasts from 4DEnVar and E3DVar analyses in these experiments show improvements over the cycling cases discussed in section 4, which include the removal of a northward track bias and the prediction of Karl’s genesis before making landfall on the Yucatan Peninsula (cf. Figs. 7 and 8). Nevertheless, these forecasts contain a slow bias for the timing of genesis that is similar to the delayed intensification found in the control simulations (Fig. 8b). The 4DEnVar-ep case yields marginal improvements over 4DEnVar, including more accurate intensity predictions in forecasts generated at earlier analysis times, but the advantages become negligible for analyses close to the genesis time. The E4DVar analyses produce a larger reduction in the intensity bias over the control, demonstrating additional skill over the other methods in forecasting Karl’s evolution. Furthermore, the E3DVar and 4DEnVar analyses produce a larger amount of variability in track forecasts between data assimilation cycles compared to the E4DVar and control simulations. The increased variability in track is partly due to the higher magnitude and variance of analysis increments generated over these cycles (see discussion at the end of this section).
To summarize differences in the four sets of analyses, we again compare time series of 1°-averaged 950- and 500-mb ζ, vortex tilt, shear, and 950–500-mb CRH near the storm center. These quantities are plotted for each assimilation cycle between 1800 UTC 12 September and 0000 UTC 15 September in Fig. 9, along with values from the control state, or prior mean used to perform the data assimilation. The time series show the low- and midlevel circulation in the 4DEnVar and E4DVar analyses being comparable to the prior mean for most cycles of the experiment. 4DEnVar-ep tends to produce a stronger low-level vortex between 13 and 14 September, which likely causes the faster intensification of Karl in these cycles compared to 4DEnVar. Likewise, E3DVar produces the weakest 950- and 500-mb circulation from 14 September onward, which is consistent with analyses from the cycling data assimilation experiments (see Figs. 3a,b). Another result that is consistent with the cycling experiments is the persistently higher values of CRH in the E4DVar analyses compared to analyses produced by 4DEnVar and E3DVar (Fig. 9e). These experiments confirm that the higher values of saturation found in the E4DVar analyses, compared to the other experiments, come mostly from the use of a tangent linear model and its adjoint to model temporal covariances. Another factor leading to differences in these results is the microphysics schemes used by the nonlinear and tangent linear models: 4DEnVar uses the nonlinear WRF Model with single-moment 6-class microphysics to approximate covariances, while E4DVar uses a simplified large-scale condensation microphysics scheme in the tangent linear and adjoint models. 4DEnVar-ep also produces conditions that are more favorable for genesis than 4DEnVar, including a decrease in vortex tilt and increase in mean ζ and CRH near the storm center at some cycles. Nevertheless, the column of air in the vicinity of the 4DEnVar-ep vortex remains consistently less saturated than the E4DVar case, which is one reason why Karl’s intensification rate in forecasts contains a slow bias that is similar to the 4DEnVar experiment at most cycles.
The time series of analysis quantities in Fig. 9 show differing magnitudes for analysis increments generated by each method during the 10 cycles examined in this section. To compare the evolution of these increments with time, Fig. 10 shows the mean squared difference between the control simulations and forecasts generated from each analysis; values are plotted for the 850-mb ζ (left column) and CRH (right column) from 0 to 72 h. The increments are separated according to wavelength (L) using high- and low-pass filters, then binned to reflect increments at the system scale (L > 150 km), intermediate scale (50 < L < 150 km), and convective scale (L < 50 km). The mean squared increments are averaged spatially over a 2000 × 2000 km2 region near the storm location, and then averaged over the 10 forecasts in this sample (solid lines in Fig. 10). Likewise, dashed lines in Fig. 10 show the standard deviation of the mean squared increments.
On average, E4DVar produces the smallest ζ analysis increments for all three sets of filtered data. At the system scale, the average squared ζ differences between E4DVar and control simulations is less than half the values found in the other data assimilation experiments at the analysis time. The E4DVar increments grow at a faster rate, and approach the same magnitude as increments from the other experiments by 36 h (Fig. 10a). A similar result occurs for the intermediate- and convective-scale ζ field (Figs. 10b,c), except the rate at which the increments reach saturation is increased for the smaller scales. The growth of CRH increments with time is consistent with the ζ increments for the intermediate and convective scales because of the relationship between moist processes and the intensification of the storm. The system-scale CRH increments, however, have negative or zero growth during the first 6 h of the 4DEnVar-ep and E3DVar forecasts, suggesting that portions of the moisture and temperature increments made by these data assimilation methods do not persist through the forecast period. Both of these experiments also exhibit higher variance in the mean squared CRH increments at the initial time, suggesting a possible overfitting of the observations. E3DVar produces larger analysis increments because the assimilated observations are binned over a longer time window, thus leading to more substantial innovations. Likewise, 4DEnVar-ep produces larger analysis increments than the standard cost function form of 4DEnVar because the static component of the background error covariance is able to grow over the assimilation window via the hybrid ensemble perturbations. The increase in background error variance in 4DEnVar-ep leads to larger analysis increments from this method. The E4DVar and 4DEnVar experiments yield system-scale analysis increments to CRH that have a similar magnitude, and exhibit a log-linear increase with time. Nevertheless, the E4DVar increments (for all wavelengths) approach a larger mean and variance by 24 h in the forecast period because of the faster intensification rate of the storm in these simulations.
This study uses three ensemble-variational hybrid data assimilation methods to generate analyses and forecasts for Hurricane Karl, an Atlantic storm targeted during the PREDICT and GRIP field campaigns. Three variational schemes are coupled with an EnKF data assimilation system, which generates posterior ensemble perturbations that are added to the hybrid variational analysis during each data assimilation cycle. The methods include hybrid E3DVar and E4DVar, as well as a hybrid 4DEnVar scheme that combines 3DVar-FGAT with temporal covariances estimated from an ensemble of model trajectories. Hybrid E4DVar and 4DEnVar apply two different strategies for assimilating observations over a time window: E4DVar uses linearized models to propagate implicitly in time a climate-based covariance and a localized ensemble covariance from the beginning of the assimilation window, while 4DEnVar uses a static climate-based covariance and localized covariances from a previously run ensemble forecast through the assimilation window. E4DVar has the advantage of propagating a full-rank hybrid background error covariance in time, while 4DEnVar benefits from using the full nonlinear model to estimate temporal covariances and is more computationally efficient to run on large parallel computing platforms. To the best of our knowledge, this study presents the first comparison of these two ensemble 4D data assimilation techniques for a tropical cyclone application.
The data assimilation systems are cycled over a weeklong period to assimilate routinely collected observations and dropsonde measurements from the PREDICT field campaign. The cycling period spans 4 days leading up to Karl’s genesis, as well as its initial intensification, and deterministic forecasts from each hybrid analysis are run over the entire life cycle of the storm. The meso-α-scale vortex in the E3DVar, E4DVar and 4DEnVar analyses contains qualitatively similar kinematic structure throughout the cycling period, including comparable low- and midlevel circulation strength and tilt. The two 4D data assimilation systems, however, produce values of column saturation that are more favorable for genesis. Elevated relative humidity values found near the developing storm in E4DVar and 4DEnVar analyses coincide with a faster strengthening of the mid- and low-level mesoscale vortex in the analyses. The circulation region near the pregenesis disturbance approaches saturation faster in the E4DVar experiment than in the 4DEnVar experiment, which coincides with earlier genesis times in forecasts. This factor leads to more accurate predictions for Karl’s formation before making its first landfall over the Yucatan Peninsula.
To remove the effects of errors that accumulate during successive data assimilation cycles, the three methods are compared using the same background ensemble in the variational analyses. These experiments are performed over a period that spans the pregenesis phase of Karl’s development and initial intensification, revealing systematic differences in the mesoscale structure of the storm. Analyses produced from these experiments show a stronger circulation in the low-level (950 mb) vortex in the E3DVar and 4DEnVar cases, compared to E4DVar. Nevertheless, E4DVar analyses contain the highest values of column relative humidity in a deep layer near the disturbance center, which contributes to faster rates of intensification in forecasts. E4DVar also produces systematically smaller analysis increments than E3DVar and 4DEnVar, which grow faster in time during the first 0–12 h of the forecast period. This result suggests E4DVar targets different dynamical instabilities during 4D data assimilation than 4DEnVar, thus reflecting differences in how covariance is modeled by the two methods.
The current study highlights several differences between hybrid E4DVar and 4DEnVar for a weather application where mesoscale processes play an important role in the system’s development. Additional data assimilation experiments and a more idealized model setup are necessary to fully explore some of the findings in this study. In particular, the modeling of covariances related to moist processes appears to be a major difference between E4DVar and 4DEnVar. We hypothesize that E4DVar may provide better estimates of temperature and moisture covariances because of its superior localization strategy. E4DVar may also benefit from using a simple linearized microphysics scheme, which improves the conditioning of the minimization problem. This study motivates further testing of E4DVar and 4DEnVar using a diverse set of case studies to gain a more comprehensive understanding of the advantages and deficiencies of these methods for mesoscale data assimilation.
This work was supported in part the Office of Naval Research Grant N000140910526 and the National Science Foundation Grant ATM-0840651. The first author is also supported by a National Center for Atmospheric Research Advanced Study Program fellowship. Computing for this study was performed at the Texas Advanced Computing Center. The authors thank Juanzhen Sun, Altug Aksoy, and three anonymous reviewers for their helpful comments on the manuscript.
The National Center for Atmospheric Research is sponsored by the National Science Foundation.
See section 3 of Poterjoy and Zhang (2015) for a discussion on the different localization strategies used by 4DEnVar and E4DVar.
In practice, the cost function is preconditioned by replacing with , where is a square root of the background error covariance and is the new control variable. For simplicity, this substitution is not shown in the manuscript.
We enforce this constraint for simplicity only. Bishop and Satterfield (2013) show analytically that the optimal weighting coefficients do not have to follow this constraint.
Another option is to use increments from the middle of the assimilation window. Section 4f of Poterjoy and Zhang (2015) provides an explanation as to why the approach taken in the current study is more accurate than the alternative method when observations are spaced equally in time through the assimilation window.
We use an 800-km radius around the best-track center for this analysis because this region includes nearly all the PREDICT and GRIP dropsondes at each verification time.