Correcting Fast-Mode Pressure Errors in Storm-Scale Ensemble Kalman Filter Analyses

A typical storm-scale ensemble Kalman filter (EnKF) analysis/forecast system is shown to introduce imbalances into the ensemble posteriors that generate acoustic waves in subsequent integrations. When the EnKF is used to research storm-scale dynamics, the resulting spurious pressure oscillations are large enough to impact investigation of processes driven by nonhydrostatic pressure gradient forces. Fortunately, thermodynamic retrieval techniques traditionally applied to dual-Doppler wind analyses can be adapted to diagnose the balanced portion of an EnKF pressure analysis, thereby eliminating the fast-mode pressure oscillations. The efficacy of this approach is demonstrated using a high-resolution supercell thunderstorm simulation as well as EnKF analyses of a simulated and a real supercell.


Introduction
The EnKF [1] has become a popular and valuable tool for storm-scale research [2][3][4][5][6][7][8][9][10][11][12].Particularly when dual-Doppler radar data are available, EnKF data assimilation can provide reliable analyses of wind and, to a lesser degree, temperature and microphysical variables in convective storms.EnKF analyses of pressure, on the other hand, are subject to severe errors, at least with some compressible model configurations (the first tests of the EnKF with a compressible model were performed by Tong and Xue [5]).This problem is illustrated in Figure 1 using output from the National Severe Storms Laboratory Collaborative Model for Multiscale Atmospheric Simulation (NCOMMAS; [13,14]) ensemble square root filter [15].Similar behavior occurs using the Data Assimilation Research Testbed [16] EnKF with the Advanced Research Weather Research and Forecasting (WRF-ARW; [17]) model (James Marquis and Thomas Jones, personal communication 2013).The pressure analysis errors severely impede investigation of critical storm processes that are, in part, driven by dynamic pressure gradient forces, including supercell occlusion downdrafts [18], lifting of negatively buoyant air [19], supercell propagation [20], the descending rear inflow and ascending front-to-rear flow in mesoscale convective systems [21], and possibly descending reflectivity cores [22].
The pressure oscillations shown in Figure 1 are associated with acoustic waves generated within the data assimilation region.The waves are presumably excited as each ensemble member adjusts to an updated initial condition that is dynamically inconsistent with the model (i.e., unbalanced).This hypothesis is supported by two observations.First, the acoustic waves occur whether or not pressure is updated during the data assimilation and therefore cannot be attributed to erroneous ensemble covariances between the pressure and other variables (though the latter could conceivably exacerbate the problem in cases where pressure is updated).Second, spurious waves are not evident in perfect-model EnKF observing system simulation experiments (OSSEs) with the NCOMMAS (not shown), indicating that the waves arise only when analysis increments are substantially unbalanced.The generation of high-amplitude fast modes due to unbalanced initial conditions is a long-recognized problem in numerical weather prediction, and many approaches have been used to improve dynamical balance during the data assimilation  procedure [24,25].Given the small influence of the pressure field on the remaining state variables in certain compressible cloud models on the relevant spatiotemporal scales, errors due to the acoustic waves are largely confined to the pressure field in storm-scale EnKF analyses.This makes it possible to retrieve the portion of the pressure field that is in balance with the remaining model fields (hereafter, the "balanced" pressure).This obviates the need to mitigate the acoustic waves during the data assimilation procedure, at least when the analysis is not being used to initialize a numerical forecast.Instead, the balanced pressure can be retrieved after the data assimilation is complete.Several methods are available for diagnosing perturbation (from the hydrostatic base state) pressure, typically in its nondimensional Exner function form  (the prime symbol is omitted herein to simplify the notation), from the equations of motion.Early methods satisfied the horizontal equations of motion on individual horizontal planes [26,27].A major limitation of that approach is that the analyzed  is offset from the true  by a vertically varying constant, precluding unique solution of / unless independent pressure measurements are available at each analysis level.This problem is avoided when all three equations of motion are satisfied [28][29][30], in which case the retrieved  is offset by a volume-wide (rather than vertically varying) constant, thereby permitting the impact of vertical pressure gradients to be quantitatively considered.The tradeoff is that errors in the analyzed local derivative and buoyancy terms in the vertical momentum equation typically result in / and / analyses that are inferior to those obtained using twodimensional retrieval.As will be shown, however, there is a very simple procedure for obtaining the advantages of both the 2D and 3D approaches.This combined method permits useful pressure retrievals to be obtained from acoustic wavecontaminated EnKF analyses.

Pressure Retrieval Method
We adopt a variational framework for our pressure retrieval scheme.The momentum equation constraints use the Klemp and Wilhelmson [31] formulation of the equations of motion except with the Coriolis term (which is negligible) omitted.
A horizontal smoothness constraint is imposed to filter noise.The cost function we seek to minimize can therefore be expressed as where , , and  are the model grid indices; , , and  are the sums of the local time derivative, advection, and turbulent mixing terms for the -, V-, and -equations, respectively;  is the gravitational acceleration; and  0 and   are the horizontally homogeneous base state and perturbation potential temperature, respectively.The subscripted 's represent the weighting coefficients for each respective constraint, computed similarly to the coefficients in [30]: where Δ, Δ, and Δ are the analysis grid spacings and  0 is the cutoff wavenumber, for which the theoretical filter response is 0.5.The present study uses a 3Δ cutoff wavelength.
In the first step of the pressure retrieval scheme, a 2D retrieval is performed (i.e.,  3 is set to zero) over each analysis level.The resulting analysis,  2D , is stored, and then a 3D retrieval is performed, yielding  3D .The vertically varying constant () by which  2D is offset from the balanced ,  bal , is then estimated by computing the mean difference between  3D and  2D for each .To see why () can be estimated in this way, consider that the 2D and 3D analyses can be written as where  error 2D and  error 3D are the retrieval errors, apart from the vertically varying constant.Assuming that the 2D and 3D retrievals are approximately unbiased relative to each other (apart from the vertically varying constant), averaging the differences between the retrievals over each horizontal analysis plane (symbolized by brackets) yields The final analysis  final is then obtained by subtracting the () estimate from  2D : The  final is still subject to a volume-wide constant, which can be estimated as the volume-mean difference between  final and the  from the simulation or EnKF analysis to which the retrieval method is applied.This estimate is reasonable to the extent that the errors in the input  are unbiased.Errors in the estimate of this constant will not impact most applications since it is typically only the spatial derivatives of  that are sought.The estimated volume-wide constant has been subtracted from the  final in the experiments below.
To further facilitate interpretation,  has been converted to dimensional perturbation pressure .

Verification of Retrieval
Procedure.An NCOMMAS supercell simulation was used to test the robustness of the retrieval procedure and the integrity of the computer code.
The simulation was performed on a stationary 102.6 km × 102.6 km × 20.4 km domain with 600 m grid spacing in all three dimensions and a large (small) time step of 4 (2/3) s.A fully dual-moment version of the Ziegler [32] microphysical parameterization (Ziegler Variable Density or ZVD) scheme [33] was used.The remaining model settings were identical to those of Potvin et al. [23].The supercell in the present simulation appears reasonably realistic, and its evolution qualitatively resembles that of Potvin et al. [23].Pressure retrievals were performed over the entire simulation domain.The local derivatives of , V, and  were computed using centered finite differences with default Δ = 30 s (e.g., / = (( + 30) − ( − 30))/60).In general, spatial derivatives were computed using centered differencing on the Arakawa A (unstaggered) grid, necessitating averaging of the model , V, and  from the Arakawa C (staggered) grid.At the model boundaries, the  derivatives in  1 −  3 were computed using one-sided differences.In preliminary retrievals, setting the turbulence terms in the equations of motion to zero generally had negligible or mildly positive impact, presumably due to discretization errors being of similar magnitude to the turbulence terms themselves, which are typically much smaller than the remaining momentum equation terms.The turbulence terms were consequently omitted in the experiments shown.
Visual comparisons of the model and retrieved  final ,  final /, and  final / reveal high fidelity in the retrieval technique (a representative example is shown in Figure 2), as do vertical profiles of root mean square error (RMSE) and relative (i.e., as a percentage of RMS  true ) RMSE (RRMSE) computed within the storm (Figure 3).(Both the model and retrieved pressure gradients were computed using centered differences on the Arakawa A grid.) Calculations of the horizontal domain-wide  bias at each level indicate that the 3D retrieval largely eliminates the vertically varying constant present in the 2D retrieval (Figure 4).To assess the sensitivity of the technique to hydrometeor mixing ratio errors (which are often large in EnKF analyses due to gross imperfections in current microphysical parameterization schemes), the retrieval was repeated with the hydrometeor drag term (in the vertical equation of motion) omitted.While ignoring water loading substantially degraded  3D within the storm (Figure 5),  3D was barely impacted in precipitation-free regions (not shown), resulting in relatively small domain-wide  3D biases at each level and, thus, only minor error increases in the () estimates (Figure 4).0.05 hPa km −1 ) by the degraded () was practically negligible.The relative insensitivity of the proposed retrieval method to hydrometeor errors is a major advantage over existing 3D approaches.The retrievals were also relatively insensitive to temporal discretization errors in the local wind derivatives.Computing the latter using Δ = 120 s (rather than 30 s) to match typical storm-scale EnKF analysis cycle periods substantially increased the RMSE  final / and  final / (Figure 3) but did not seriously qualitatively degrade the retrieval (Figure 6).While temporal discretization errors will be much larger in cases of rapidly moving storms (the present storm travels eastward at ∼10 m s −1 ), such errors may be substantially reduced using advection correction methods [34][35][36].

Retrievals from a Simulated EnKF Mean Analysis.
Having verified the robustness of the pressure retrieval scheme, we then applied the technique to mean posteriors from an EnKF OSSE conducted by Potvin and Wicker [11] (pictured in Figure 1).In that experiment (labeled "2-LFO" in [11]), mobile dual-Doppler velocity pseudoobservations were generated from the "truth" simulation of Potvin et al. [23], perturbed with random errors, and assimilated at two-minute intervals.Model error was introduced by using coarser ensemble grid spacing (600 m versus 200 m) and different microphysical parameterization (the Gilmore et al. [37] version of the Lin et al. [38] scheme versus the ZVD scheme) from those in the "truth" simulation.The use of an imperfect model in the EnKF resulted in dynamically unbalanced analysis increments that instigated acoustic waves during the data assimilation.The ensemble was integrated for 30 min prior to the first analysis update to develop physically realistic covariance structures.
Local wind derivatives for the pressure retrieval scheme were computed using mean posteriors separated by 4 min (i.e., Δ = 2 min).This interval is similar to the volume scan periods for Weather Surveillance Radar-1988 Doppler (WSR-88D) convection sampling patterns.All fields required for the momentum equation constraints were obtained from the EnKF posteriors (i.e., treated as known quantities).The pressure retrievals were performed offline from the EnKF analysis; that is, the retrieved posterior pressure field at each time was not used to initialize the subsequent forecast cycle in the data assimilation experiment.Due to the very weak impact of the pressure field on the remaining state variables in compressible models (the property exploited by the proposed method to remove the fast-mode pressure oscillations), rebalancing the pressure field at each analysis cycle should have generally negligible impact on the remaining model fields throughout the data assimilation.We have verified this with separate, real data experiments (not shown), in which the differences due to the pressure rebalancing produced only tiny differences in the final ensemble analysis and in subsequent 30 min forecasts.Given the insensitivity of the pressure retrievals in Section 3.1 to discretization and hydrometeor mixing ratio errors, it was expected that  final obtained from the EnKF analyses would well capture the balanced portion of the EnKF  ( bal ) and thereby improve upon the EnKF  and generally resemble  true (keeping in mind that errors in the other EnKF fields will create some discrepancies between the EnKF  bal and  true ).This was indeed the case throughout the assimilation period, except at higher altitudes.After 20 min of data assimilation ( = 50 min), the EnKF wind analyses are sufficiently accurate that  final reasonably resembles  true at lower and middle levels, much more so than does the EnKF  (Figures 7, 8, and 9).The improvement of  final over the EnKF , however, decreases with height (Figure 9), presumably due in large part to the shallower fast-mode pressure oscillations aloft (not shown).Within the top half of the storm,  final is generally mildly inferior to the EnKF , indicating that retrieval errors dominate improvements from removing the (small) fast-mode pressure errors.This suggests that application of the proposed pressure retrieval method be restricted to altitudes where spurious pressure oscillations are large.While the retrieval results from this OSSE are likely somewhat optimistic despite our efforts to mitigate the "identical twin" problem, larger EnKF analysis errors in real cases should not substantially impede the technique's ability to recover the EnKF  bal and, thus, to remove the fastmode pressure errors (which is the objective of the proposed method).Furthermore, the previous results, combined with the resilience of dual-mobile-Doppler EnKF wind analyses to microphysical and background wind errors [11,39], suggest that pressure retrievals obtained from high-quality radar datasets such as those collected during the second Verification of the Origins of Rotation in Tornadoes Experiment (VORTEX2; [40]) may be useful for investigating pressuregradient-driven storm processes.Rigorous testing of this hypothesis, unfortunately, is currently precluded by the lack of dense 3D pressure observations within storms.

Retrievals from a Real Data EnKF Mean Analysis.
While our ability to verify the real-world performance of the proposed technique is quite limited, qualitative evaluation of pressure retrievals from EnKF analyses of the May 29-30, 2004, Geary, OK, USA, tornadic supercell supports cautious optimism.The EnKF analyses were obtained by assimilating quality-controlled Doppler velocity data from two Shared Mobile Atmospheric Research and Teaching (SMART; [41]) radars using the NCOMMAS EnKF with the ZVD microphysics scheme (see Potvin et al. [39] for further details).Pressure retrievals were performed in the same manner as with the simulated EnKF analyses in Section 3.2.
A representative comparison of the  final and EnKF , valid after 44 min of data assimilation (0037 UTC 30 May), is shown in Figure 10.As in the EnKF OSSE (Section 3.2), highamplitude pressure oscillations generated during the data assimilation are not evident in the retrieved pressure fields.Furthermore, the  final comports with present understanding of the pressure distribution within supercells, which suggests retrieval errors are not unduly large.For example, both the  final and EnKF  exhibit an inflow low and a deep pressure minimum associated with the strong mesocyclonic rotation (Figures 10(a) and 10(b)), as well as upward-directed perturbation pressure gradient force along the gust front (Figures 10(c) and 10(d)).The  final , but not the EnKF , persistently indicates a region of substantial downward-directed perturbation pressure gradient force near the analyzed occlusion downdraft (Figures 10(c)-10(f)).This is consistent with current understanding of occlusion downdraft formation [42], suggesting that the EnKF pressure analysis near the occlusion downdraft is contaminated with fast-mode errors that the retrieval successfully corrects.

Conclusion
As in other data assimilation frameworks, storm-scale EnKF radar data assimilation can introduce dynamical imbalances that induce severe pressure errors in analyses, inhibiting investigation of important storm processes.Retrieval procedures that have traditionally been used to retrieve the pressure and buoyancy fields from dual-Doppler wind analyses can be used to rebalance the pressure field in such cases.Using a combination of the 2D and 3D pressure retrieval approaches that mitigates the deficiencies of both, we have demonstrated that useful analyses of the perturbation pressure field and its derivatives can be obtained from EnKF analyses despite model and discretization errors.The described procedure could be modified to separately retrieve, for example, the linear and nonlinear components of the dynamic perturbation pressure [43].The method could also be extended to other nonhydrostatic equation sets including that used by the WRF-ARW model.
It should be noted that some of the existing methods for suppressing acoustic and gravity waves in model initial conditions could be adapted to removing pressure oscillations in EnKF analyses.For example, during data assimilation, digital filter initialization (DFI; [44]) could be applied to the analysis increments at each cycle prior to integrating the ensemble forward to the next time.All such approaches, however, would increase computational cost, and it is not clear whether storm-scale analyses would improve given the model approximations used in these methods (e.g., the adiabatic backward model integration in DFI).Alternatively, a diagnostic pressure equation (e.g., [45]) could be solved (either exactly, or weakly as in the proposed method), thereby avoiding temporal discretization errors, but at the cost of increased spatial discretization errors.While we recommend evaluation of the relative efficacy of alternative methods for recovering the balanced pressure component, we also emphasize that the proposed method provides an existing, effective, simple, and low-cost way to remove pressure oscillations from storm-scale EnKF analyses and does not require data assimilation experiments to be rerun.

Figure 1 :
Figure1: True (left column) and EnKF mean posterior (right column)   (shading; hPa) and radar reflectivity factor (contoured at 10, 30, and 50 dBZ) at  = 0.9 km.Fields are valid after (top row) zero, (middle row) one, and (bottom row) fifteen 2 min forecast cycles.The true   were filtered and averaged as in Potvin et al.[23] to permit more direct comparisons with the (coarser) EnKF   .

Figure 6 :
Figure 6: As in Figure 2 except for Δ = 120 s in the pressure retrieval.

Figure 9 :Figure 10 :
Figure 9: As in Figure 3 but for  final retrieved from mean EnKF posterior (solid) and mean EnKF posterior  (dashed).