Wind extraction potential from ensemble Kalman filter assimilation of stratospheric ozone using a global shallow water model

The feasibility of extracting wind information from stratospheric ozone observations is tested using ensemble Kalman filter (EnKF) data assimilation (DA) and a global shallow water model that includes advection of an ozone-like tracer. Simulated observations are created from a truth run (TR) that resembles the Northern Hemisphere winter stratosphere with a polar vortex disturbed by planetaryscale wave forcing. Ozone observations mimic sampling of a polar-orbiting satellite, while geopotential height observations are randomly placed in space and time. EnKF experiments are performed assimilating ozone, height, or both, over a 10-day period. The DA is also implemented using two different pairs of flow variables: zonal and meridional wind (EnKF-uv) and stream function and velocity potential (EnKF-ψχ). Each experiment is tuned for optimal localization length, while the ensemble spread is adaptively inflated using the TR. The experiments are evaluated using the maximum wind extraction potential (WEP). Ozone only assimilation improves winds (WEP = 46 % for EnKF-uv, and 58 % for EnKF-ψχ), but suffers from spurious gravity wave generation. Application of nonlinear normal mode initialization (NMI) greatly reduces the unwanted imbalance and increases the WEP for EnKF-uv (84 %) and EnKF-ψχ (81 %). Assimilation of only height observations also improved the winds (WEP = 60 % for EnKF-uv, and 69 % for EnKF-ψχ), with much less imbalance compared to the ozone experiment. The assimilation of both height and ozone performed the best, with WEP increasing to ∼ 87 % (∼ 90 % with NMI) for both EnKF-uv and EnKF-ψχ , demonstrating that wind extraction from ozone assimilation can be beneficial even in a data-rich environment. Ozone assimilation particularly improves the tropical winds, which are not well constrained by height observations due to lack of geostrophy.


Introduction
A key missing component of the global observing system (GOS) is measurement of the three-dimensional global wind (World Meteorological Organization, 2000).Upper air wind observations from radiosondes, pilot reports, and cloud and water-vapor feature tracking leave large gaps, particularly in the tropics, Southern Ocean, and in most of the stratosphere and mesosphere.Spaceborne Doppler wind lidar (DWL) has been proposed as the potential "missing link" in the GOS (Baker et al., 2014).When placed in low earth orbit, DWL can provide daily global wind profiles throughout the troposphere and lower stratosphere (National Research Council, 2007).The Atmospheric Dynamics Mission (ADM-Aeolus) (Stoffelen et al., 2005) will provide a proof of concept of this capability.However, the measurements will be limited to a single line-of-sight wind component, altitudes below ∼ 26 km, and simple along-track (as opposed to sweeping or conical) sampling.While future spaceborne DWL missions may provide improved observing capabilities, the technical challenges make this a very difficult and expensive solution to the problem of inadequate wind observations.
Another potential source of dynamical information comes from assimilation of trace gas (tracer) observations in a 4-D data assimilation system (DAS) that dynamically couples tracer and wind.The investigation of algorithms to extract wind information from tracers started with 1-D and 2-D simulations by Daley (1995Daley ( , 1996) ) and Riishøjgaard (1996).These studies showed that wind information could be extracted from tracer observations when the continuity equation was coupled to the dynamical equations via either a 4D-Var (four-dimensional variational assimilation) algorithm or an extended Kalman filter (EKF).Extensions to the full 3-D atmosphere were performed in 4D-Var experiments by Peuch et al. (2000), Semane et al. (2009), and Allen et al. (2013).These further supported the potential of tracer assimilation to benefit the winds but also highlighted limitations due to paucity of observations, insufficient data quality, and inadequate modeling of tracers in the forecast model, as well as phenomenological limitations due to geophysical variability.
Assimilation of infrared and microwave humidity channels from geostationary and polar-orbiting satellites has been shown to benefit tropospheric analyses and forecasts in the European Centre for Medium-Range Weather Forecasts (ECMWF) 4D-Var system (Andersson et al., 2007;Peubey and McNally, 2009).Peubey and McNally (2009) isolated the mechanisms whereby geostationary clear-sky radiances can impact the wind analyses in 4D-Var and showed that the dominant factor involves adjustment of the wind field to match observed humidity features (the so-called "tracer advection effect").However, attempts to assimilate stratospheric ozone using 4D-Var algorithms and the resultant dynamical coupling have previously resulted in problems in operational numerical weather prediction (NWP) (Han and McNally, 2010;Dragani and McNally, 2013).These assimilation challenges led Allen et al. (2014) to re-examine the stratospheric tracer-wind problem at a more fundamental level using 4D-Var assimilation studies with a shallow water model (SWM) coupled to the tracer continuity equation.This idealized system allowed Allen et al. (2014) to probe the limits of wind extraction from assimilation of three readily measured long-lived tracers: ozone (O 3 ), nitrous oxide (N 2 O), and water vapor.It was shown that assimilation of global hourly tracer data was sufficient to analyze the horizontal wind components to a high degree of accuracy (∼ 0.3 m s −1 random error for O 3 and N 2 O).
While 4D-Var couples tracers and dynamical variables through the tangent linear model and its adjoint, the initial background error covariance normally does not include tracer-wind correlations (these correlations develop implicitly over the assimilation window).This limitation may be overcome by using an ensemble Kalman filter (EnKF) in which the error covariance between tracer and wind is explicitly calculated by the ensemble statistics.Milewski and Bourqui (2011) assimilated ozone and temperature profiles in an EnKF system using a 3-D model at relatively low resolution (spectral triangular truncation T21).They showed that background error covariances are able to propagate information from the observed variables to wind.In particular, assimilation of either ozone or temperature observations in a polar-orbiting sampling pattern significantly improved the wind analysis.Another approach to enhancing the tracerwind interaction within 4D-Var is to blend the static covariance with a flow-dependent ensemble covariance.This hybrid 4D-Var method is becoming increasingly popular at operational NWP centers (Buehner et al., 2010;Bonavita et al., 2012;Clayton et al., 2013;Kuhl et al., 2013;Kleist and Ide, 2015).We are developing a hybrid system within the SWM framework to study tracer-wind interaction, which we plan to present in a follow-up paper.
In this paper, we take a similar approach to Milewski and Bourqui (2011), except that we use the SWM forecast model (at T42 resolution), and we assimilate ozone and height (in lieu of temperature for the SWM) observations, both separately and together, to examine whether value is added by assimilating ozone observations into a system already constrained by other observations.The SWM has been used in both 4D-Var (Courtier and Talagrand, 1990;Polavarapu et al., 2000;Jung et al., 2014) and EnKF (Kepert, 2009(Kepert, , 2011) ) experiments, since it provides a sufficiently complex system to simulate the key physical relations of the horizontal flow, including both slow balanced and fast unbalanced modes.As explained by Kepert (2009), the SWM provides a severe test for assimilation, since the weak dissipation will not remove imbalances introduced in the analysis; they will rather accumulate with time.
One of the goals of the current study is to probe the limits of ozone-wind extraction in an EnKF system.To accomplish this, it is necessary to quantify (and remove, if possible) spurious imbalance generated from noisy observations and imperfect modeling of background error covariances.A wide range of studies has been performed to examine balance in the context of 4-D data assimilation.For example, Neef et al. (2006Neef et al. ( , 2009) ) investigated balance with a loworder Lorenz-type model with the EKF and EnKF.Imbalance within SWM-DAS systems was analyzed in both 4D-Var (Courtier and Talagrand, 1990;Polavarapu et al., 2000) and EnKF (Kepert, 2009(Kepert, , 2011) ) using digital filter and nonlinear normal mode initialization techniques.Mitchell and Houtekamer (2002) considered the influence of covariance localization on balance with a 3-D dry, global, primitiveequation model.In all of these studies, imbalance was shown to be a serious issue in 4-D DAS.None of these studies was designed to examine balance in the context of the tracer assimilation problem, however.As part of this study we attempt to isolate and to the extent possible remove imbalance in order to minimize the analysis errors and determine the extent to which the wind can be constrained by ozone observations.The layout of the paper is as follows.Section 2 describes the SWM-DAS, including the forecast model, the EnKF, and the normal mode initialization procedure.Section 3 describes the experimental design and the error diagnostics.Section 4 presents the results and discussion from the three assimilation experiments, and conclusions are provided in Sect. 5.

Forecast model
The forecast model is a spectral SWM based on the vorticitydivergence formulation in Sect.2a of Ritchie (1988), with the inclusion of fourth-order semi-implicit diffusion applied to the vorticity, divergence, and geopotential.A spectral advection equation is coupled to the SWM, solving for the mixing ratio of a passive tracer as a function of time using the same fourth-order diffusion operator.We call the combined four-equation system the shallow water model with tracer (SWM T ).The system is run at triangular truncation T42, with model fields saved on the Gaussian grid (128 longitudes × 64 latitudes, for a grid resolution of ∼ 2.8 • at the Equator).The discretization uses a leap-frog time integration and a semi-implicit approximation for terms that produce gravity waves (Ritchie, 1988).To restart the model after assimilating data, a forward Euler time stepping method is applied.The global mean geopotential height H is specified to be 10 km, resulting in a gravity wave speed ( √ gH , where g is earth's gravitational acceleration) of 313 m s −1 .To avoid numerical instability due to gravity waves, a short model time step of 120 s is used for all SWM T forecasts.The diffusion coefficient is set to 5.0 × 10 15 m 4 s −1 , which provides an e-folding damping for the highest wave number of approximately 1 day.

Ensemble Kalman filter
To assimilate data into the SWM T system, we use the "perturbed observations" EnKF (Houtekamer and Mitchell, 1998;Evenson, 2003).The system solves for N ens analysis states using the Kalman filter equation for the state vector x of size N state .
where a and b superscripts indicate analysis and background, and i = 1...N ens is an index for ensemble member.d i = y i − Hx b i is the vector of innovations for member i, y i is the vector of perturbed observations, and H is the (linear) observation operator.The ensemble-based Kalman gain matrix is defined as follows: with the ensemble background error covariance calculated by Here the overbar indicates the ensemble mean and R is the observation error covariance matrix.The background state is calculated using the nonlinear SWM T forecast model M, subject to initial conditions where n is an index for model time (t n ).Note that the SWM T time step (120 s) is less than the analysis time step (20 min), such that there are 10 forecast model time steps between analyses.At each analysis time, which corresponds to the end of the 20 min background forecast, all the observations at that time are assimilated simultaneously as a single batch.The EnKF analysis equation can be solved using different combinations of state variables.In this study, we compare results using zonal wind, meridional wind, height, and ozone x = u, v, z, q (the EnKF-uv system), and results using stream function, velocity potential, height, and ozone x = ψ, χ , z, q (the EnKF-ψχ system).The latter combination was shown by Kepert (2009) to result in better balance of increments in a SWM-EnKF system with Schur product localization (discussed further below).We will test this for the SWM T system with ozone and height observations.
To avoid filter divergence, we apply a state space covariance inflation factor (Anderson, 2007) to the background ensemble before assimilating observations.The background ensemble perturbations (x b i − x b ) are multiplied by a scalar factor ξ to produce the inflated background ensemble which alters the background error covariance (Eq.3), but leaves the background ensemble mean, x b , unchanged.
The inflation factor is designed to alter the global average SPREAD to match the global root mean square error (RMSE) of either the vector wind (for EnKF-uv) or the stream function (for EnKF-ψχ ).The RMSE and SPREAD for vector wind are defined as where V represents the magnitude of the vector wind and TR indicates the truth run (described in Sect.3.1).For stream function (ψ) the RMSE and SPREAD are defined as The inflation factor is defined as either ξ V = V RMSE /V SPREAD or ξ ψ = ψ RMSE /ψ SPREAD .While the inflation factor is calculated using only the wind or stream function, it is applied to the entire state vector for each ensemble member using Eq. ( 5).The calculations of RMSE and SPREAD are not area-weighted and therefore may be somewhat biased to match the higher latitudes, since the Gaussian grid is used.This tuning takes a similar approach to the 4D-Var simulations of Allen et al. (2014) in which the background error variances were modified to match the global RMSE of the tracer and wind components.This adaptive tuning provides a flexible way to examine how the system behaves over a wide range or parameters, without needing to separately tune the inflation factor for each case.It is of course not practical in an operational setting, since the true state is unknown, but for this idealized study it works well to prevent filter divergence.
To avoid spurious long-range correlations, localization is applied to the background error covariance.We apply the element-wise (Schur product) approach (e.g., Houtekamer and Mitchell, 2001) using Eq.(4.10) of Gaspari and Cohn (1999).The localization matrix S is applied directly to the background error covariance so the gain matrix becomes To illustrate the ozone-wind interaction in the SWM T -EnKF system, Fig. 1 shows the ensemble mean analysis increments K ens d for assimilation of a single ozone observation at 120 • E longitude, 40 • N latitude for the EnKF-uv and EnKFψχ systems.For the EnKF-ψχ, we convert the increments to wind increments after the analysis step.Specification of the initial 100-member ensemble for this system will be discussed in Sect.3. The positive ensemble mean ozone innovation ( d = 0.21 parts per million by volume, ppmv) results in a positive ozone increment in the vicinity of the observation with maximum of 0.11 ppmv.Since the height correlates positively with ozone, a positive height increment also occurs (maximum of 84 m).Note that the ozone and height increments are similar for both systems, since these variables are unchanged; slight differences are due to differences in tuning of the background error covariances.The wind increments are very different, however.While both show anticyclonic circulation around the positive height increment, the winds are much stronger in the EnKF-ψχ .As explained by Kepert (2009Kepert ( , 2011)), the weakening of the winds in the EnKF-uv is due the effects of localization, which acts to decrease the local balance.As shown below, this adversely affects the system by generating spurious gravity waves.

Normal mode initialization
In general, analysis increments may project onto both slow balanced modes and fast unbalanced modes.Unless there is sufficient information in the background error covariance to distribute increments in a balanced way, the unbalanced modes will enter the system, and it may be difficult to remove these modes with limited observations (Neef et al., 2006(Neef et al., , 2009)).To quantify the imbalance in the SWM-EnKF, we use a nonlinear normal mode initialization (NMI) procedure (Machenhauer, 1977), which has been used in NWP to reduce the impact of inertia gravity waves caused by imbalance in the analysis increments (e.g., Kleist et al., 2009).While digital filter initialization (e.g., Fillion et al., 1995) is more commonly used in NWP today, NMI allows discriminating between the gravity wave and rotational wave modes, which is very useful in the SWM context.Kepert (2009) used NMI to analyze imbalance caused by localization in the SWM-EnKF framework and showed that while gravity waves can be reduced by judicious choice of balance constraints, some initialization may still be necessary in the EnKF (see also the discussion in Lorenc, 2003).
For example, the single-observation increments in Fig. 1 result in unbalanced motions in both versions of the SWM T system.Figure 2 shows the divergence anomalies due to the single-observation increments.These anomalies propagate radially outward from the observation location, as seen in these maps at 2 h intervals.Maps at later times (not shown) indicate that these oscillations propagate around the globe in ∼ 1.5 days, consistent with waves traveling at the gravity wave phase speed of this system.The EnKF-ψχ increments result in smaller divergence fields than the EnKF-uv; the maximum divergence anomaly at 1200 s for the EnKFψχ is ∼ 13 % of that caused by the EnKF-uv, consistent with less imbalance.However, initialization may still be necessary in both systems to remove this spurious gravity "noise".
To "initialize" the system (i.e., apply NMI to the analysis state vector), we first need the normal modes (NMs) of the SWM system.These were calculated using the formulation outlined in Hogan et al. (1992).The resulting NM frequencies are shown in Fig. 3 as a function of zonal wave number and mode type.Negative (positive) wave numbers indicate westward (eastward) propagating modes.These modes are separated into westward and eastward gravity wave (GW) and westward rotational wave (RW) modes.To balance the GW modes, we apply the Machenhauer (1977) condition, which reduces the time tendencies of the complex amplitudes of the modes.We apply five iterations to solve the nonlinear balance equation using a single 120 s time step for the calculation of the tendencies.We choose a linear cutoff frequency of 1.0 day −1 , which attempts to balance all traveling modes except for one eastward wave 1 GW mode (see Fig. 3a).
In this study, we apply NMI to the ensemble-meananalyzed fields only as a post-processing diagnostic to quantify the degree of imbalance.The goal is to tune the EnKF system to minimize unwanted imbalance, without having to rely on applying NMI within the DAS.One reason to avoid initialization in the EnKF cycling is that it fails to distinguish real and spurious gravity waves and can therefore potentially move the system away from the truth.Another reason is that running NMI in the EnKF would involve initializing each ensemble member separately, since different modes may be excited in each member due to the perturbed observations, which adds significantly to the computational expense.In principle, if the unbalanced modes do not interact much with the balanced components of the flow, then it should not matter whether the balancing is done before or after the assimilation.Williamson and Temperton (1981), using a multilevel global grid-point model, showed that forecasts made  with initialized data produced virtually identical results to forecasts with uninitialized data followed by initialization.This suggests that the high-frequency GWs do not interact much with the low-frequency RWs, but rather can be largely considered "noise" in the system that can, in principle, be filtered out.To test whether this is true for the system run here, we compared results using NMI cycling and NMI post-processing for the optimal runs of the three experiments examined in Sect. 4. Differences in wind extraction potential (defined in Sect.3) were ∼ 1 % or less for all runs except for height only assimilation with the EnKF-uv system, which showed an improvement of 5 % for NMI cycling over NMI post-processing.Assimilation of height observations is likely more sensitive to GW noise, which impacts the height di- rectly, while the tracer is only indirectly impacted via the divergent component of the wind, which is small compared to the rotational wind.
To illustrate the influence of NMI post-processing, Fig. 4 shows the true divergence along with the analyzed divergence with and without NMI for a sample field 2 days into an ozone assimilation run.Whereas the true divergence is rather smooth, the uninitialized divergence shows considerable noise.After applying NMI, the analyzed divergence looks much more like the truth, indicating that the noise was due largely to spurious unbalanced modes.We note that this rather large improvement from application of NMI is partly illustrating sensitivities in the SWM.Whereas in a full NWP system physical and radiative processes may dampen the gravity waves, in the SWM with weak diffusion the waves can remain in the system for a long time, unless assimilated data are at sufficient sampling frequency and precision to resolve the waves.

Truth run
The TR is designed to simulate Northern Hemisphere (NH) winter conditions in the middle stratosphere (the same TR was used in the 4D-Var tracer assimilation study by Allen et al., 2014).The initial conditions for the SWM T include zero meridional wind and a zonally symmetric zonal wind that varies with latitude (φ) as u (φ) = u max sin (2φ) in the NH and is zero in the Southern Hemisphere (SH), with u max = 60 m s −1 .The geopotential height is specified using the gradient wind balance with a global mean height of 10 km.The initial ozone is calculated using Aura Microwave Limb Sounder (MLS) ozone data (Waters et al., 1999;Livesey et al., 2011).The data are selected for a period with weak planetary wave activity (1-15 March 2011) and are interpolated to the 850 K isentropic level (approximately 32 km altitude or 10 hPa), representative of middle stratosphere conditions.The zonal mean and time mean mixing ratio as a function of latitude was calculated for this period and interpolated to the Gaussian grid.The ozone is treated as passive (i.e., no chemical source/sink) and there is no radiative interaction between ozone and dynamics.Note that in the TR there are no "restarts", since it is a continuous free-running SWM T forecast.Therefore, an assimilation cycling run, which restarts with a forward Euler step after each 20 min analysis cycle, would produce a slightly different result from the freerunning forecast (which uses leap-frog time integration for all steps after the initial Euler step), even if no data were assimilated.We could, in principle, stop and restart the TR with a forward Euler step at the regular analysis time intervals, as was done in Kepert (2009); however, test runs performed with and without restarts in the TR resulted in negligible differences.
To create a realistic scenario of the NH winter stratosphere, the TR is forced by the bottom topography being raised and lowered to simulate planetary-scale waves (as in Norton, 1994).A mountain with a height of 1250 m is created with a 20-day cycle (4 days ramping up, 12 days constant, and 4 days ramping down).The mountain is a zonal wave 1 feature that peaks at 45 • N. The topography is turned off after 20 days.Since the assimilation period corresponds to days 20-30 of this TR, there is no surface topography during the assimilation.NH maps of the ozone and height fields for the TR are provided in Fig. 5a-f.On day 20, a strong anticyclone (indicated by H) is present near 180 • longitude, resembling an "Aleutian High", with elevated ozone values.The polar vortex (indicated by L), identified by low ozone, is displaced off the pole into a comma shape.Over the next 10 days, the "Aleutian High" diminishes in strength and the vortex moves over the pole.Strong ozone advection occurs throughout this period.For example, a long tongue of lower ozone mixing ratio forms around a secondary anticyclone centered near 60 • E longitude on day 28.This dynamical scenario produced by topographic wave forcing in the SMW T provides a realistic representation of the final stages of a stratospheric minor warming (Limpasuvan et al., 2004).In the SH (Fig. 5g-l), a strong anticyclone is centered just off the pole on day 20.This anticyclone (H) propagates westward around the pole, making one cycle over this 10-day period.A weaker cyclone (L) also propagates westward around the pole opposite to the anticyclone.The ozone is advected along with these features, with relatively high (low) ozone in the anticyclone (cyclone).The westward flow in the SH is consistent with the easterly summer flow in the middle stratosphere (Andrews et al., 1987).Additional maps of potential vorticity and ozone for this TR are provided in Fig. 2 of Allen et al. (2014).

Observations
Observations are simulated by sampling the TR ozone and height fields using a bi-linear interpolation in latitude and longitude.Gaussian random error is then added with a specified standard deviation (SD).The error SD for ozone was set to 0.08 ppmv, which is 1 % of the initial global mean, while the height error SD was set to 50 m.The height error can be approximately related to stratospheric temperature error by using a climatological estimate of the equator-to-pole gradient of temperature with respect to geopotential height in the NH winter middle stratosphere of ∼ 5 K km −1 .Using this conversion factor, a 50 m error corresponds to ∼ 0.25 K.Both the ozone and height errors are smaller than those of any current operational instrument.The goal here is not to evaluate an actual observing system, but to demonstrate ozone-wind extraction in an idealized system.The observation errors are assumed to be uncorrelated, so the observation error covariance R is diagonal with elements given by the square of the error SDs.
Two sampling methods are performed (see Fig. 6).For ozone, the observation locations are taken from real ozone observations from the Aura MLS polar-orbiting satellite  (sampling frequency of ∼ 3450 observations per day).For height, pseudo-random sampling in space and time is performed to approximate the global coverage provided by microwave and infrared radiance sensors.For each successive height observation, the sampling occurs at one of 3840 latitude/longitude points on an icosahedral equal-area grid.This allows the observations to not be too clumped together and provides a way to scale upward to a global equal area grid sampling, as was used in Allen et al. (2013).We choose the average data frequency for the height observations to be the same as the MLS sampling frequency (∼ 3450 day −1 ).In time, the height observations occur randomly over 10 days.For both observation types, the observation time is assigned to the nearest 20 min interval (0, 20, or 40 min).Since the analysis is performed sequentially every 20 min, time interpolation is not necessary for observations.

Assimilation experiments
We use 100 ensemble members for all experiments in this paper.The initial ensemble perturbations are generated by sampling the TR fields at 6 h intervals (starting day 21, 0 h) and then removing the ensemble mean.The assimilation experiments begin 20 days into the TR (day 20, 0 h) with the initial ensemble defined as the ensemble perturbations added to the TR field that is offset 6 h from the initial time (i.e., day 20, 6 h), which serves as the initial ensemble mean and defines the initial analysis.This initial 6 h offset, or mismatch, between the TR and the initial background fields is the source of the initial background error.In Sect. 4 we present results from three different experiments: (1) ozone only, (2) height only, and (3) ozone and height.For each experiment, assimilation runs were done for both EnKF-uv and EnKF-ψχ using localization lengths from 1000 to 8000 km, in 500 km increments (note that ozone only experiments failed to converge at the maximum localization length of 8000 km).Tuning the length separately for each experiment and EnKF is necessary, since the DA responds differently depending on the field(s) observed and the analysis variables used.For each experiment we use the same localization length for all state variables.Further optimization may occur by applying different localization functions to different variables, but this is beyond the scope of this first study on tracer-wind interaction using the EnKF.Since inflation is automatically adjusted in a self-consistent manner with the TR, it does not require tuning.Post-processing with NMI was also performed for each run.We note here that the same forecast model is used for the TR and for the assimilation experiments (i.e., "identical twin" experiments), making results overly optimistic.

Error metrics
To diagnose the results, several error metrics are examined, including the global RMSE (area-weighted) of the u, v, z, and q, along with the wind extraction potential (WEP).Allen et al. ( 2014) defined WEP as a normalized diagnostic of the impact of tracer assimilation on the dynamics.The WEP is determined by first calculating the analyzed RMSE of the vector wind as a function of latitude and time (t): where is the RMSE of the zonal wind calculated around a latitude circle containing N lon longitude (λ) grid points and TR refers to the truth run (the RMSE of the meridional wind, v RMSE , is calculated similarly).Here i = 1... N lon is an index for longitude and j = 1... N lat is an index for latitude, both on the Gaussian grid.The latitude dependence is shown explicitly here, since we will examine errors as a function of latitude in Sect. 4. The percentage difference in vector wind error relative to the initial error is then calculated, and WEP is defined as the area-weighted global average of this quantity, calculated using where the summation is over all N lat latitude grid points.
A WEP value of 100 % indicates the analysis equals the truth (i.e., V RMSE φ j , t = 0 at all latitudes).Although WEP is relative to the initial error, and therefore will vary from one experiment design to another, it provides a useful normalized number for quantitative comparison between runs using the same initial error.In this paper, all experiments start with the same initial vector wind error with a global mean value, V RMSE (0) = 4.55 m s −1 , so WEP can be compared directly among all runs (tilde refers to the global areaweighted mean).As a rule of thumb, an approximate conversion from WEP to wind component error can be derived by assuming RMS (root mean square) wind errors do not vary with latitude and assuming zonal and meridional wind errors are equal.This results in the following approximation: The initial global mean zonal wind RMSE, u RMSE (0), is ∼ 3.3 m s −1 , so WEP values of 50, 60, 70, 80, and 90 correspond to approximate wind component errors of 1.65, 1.30, 1.00, 0.66, and 0.33 m s −1 , respectively.Experiment errors are generally presented as the "final" error of the 10-day simulations.To reduce random noise, the final errors are calculated as the average values over the last 24 h of each simulation.To estimate the statistical uncertainty in the final errors, Experiment 3 was repeated 10 times with different random observation perturbations, with a localization of 3500 km.The SD of the final values was ∼ 0.5 % for WEP, ∼ 0.02 m s −1 for the wind components, ∼ 0.4 m for height, and ∼ 0.002 ppmv for ozone.The results in Table 1 are presented with the number of significant digits that reflect the uncertainties determined from this test.
The final diagnostic is designed to measure the amount of gravity wave "noise" in the system, also called "imbalance".Imbalance is defined here using the uninitialized height z uninit and initialized height z init (i.e., after NMI has been applied).
As with WEP, we first calculate the RMS difference between these two fields as a function of latitude, and then we calculate the area-weighted global mean, We note here that the TR used in this paper contains negligible gravity wave amplitudes at frequencies higher than 1.0 day −1 .The imbalance calculated by applying NMI to the TR is less than 1 m.So any imbalance greater than ∼ 1 m is due to spurious GW generation.

Experiment 1: ozone only
In this section, we examine the performance of the SWM T -EnKF system when ozone data are assimilated alone.Figure 7 shows time series of error diagnostics for the two "optimal" runs from Experiment 1 (tuning of the covariance localization to determine the optimal run is described later in this section).Results are presented both for uninitialized (solid lines) and initialized (NMI, dotted lines) output.The uninitialized WEP steadily increases before leveling off at final values of ∼ 46 % for EnKF-uv and ∼ 58 % for EnKF-ψχ , with corresponding wind component errors of ∼ 1.6 m s −1 and ∼ 1.3 m s −1 .Most of the improvement occurs during the first 5 days.After applying NMI to these runs, the initialized WEP increases significantly for both systems, indicating that imbalance is limiting the wind improvement and error reduction.
The uninitialized height error for EnKF-uv levels out at ∼ 61 m while for EnKF-ψχ the uninitialized height error reaches ∼ 43 m.Much of the height error can be attributed to GWs generated in the system.Figure 7e shows that the imbalance starts near zero, but increases as GWs are introduced into the system.For EnKF-uv, the imbalance rises rapidly over the first 2 days until it nearly matches the uninitialized height error.After this time, the further growth of GW is likely restrained by the weak dissipation in the SWM T system.For EnKF-ψχ, the imbalance grows more slowly but is still close to the uninitialized height error at the end.Because the uninitialized height error and imbalance are still increasing at the end of 10 days, this suggests that the GWs have not saturated for EnKF-ψχ .Application of NMI results in dramatically reduced height errors for both systems.There is an ∼ 1 day oscillation in the initialized height errors, which is largely due to the one eastward-traveling GW mode that is not initialized, which is present in the TR.The decreasing amplitude of the oscillation with time suggests that this mode is being resolved by the system through the ozone-wind extraction.
The ozone errors (Fig. 7f) show a sharp decrease over the first day, followed by a gradual decline.Although the errors in the dynamical variables have leveled out by day 10, the ozone errors appear to still be declining.Both runs show final ozone errors smaller than the observation error of 0.08 ppmv.As a global consistency check of the EnKF solution, we also calculated where d is the ensemble mean innovation and N obs is the number of observations.For a well-tuned system, this "chisquared" diagnostic should equal 1 (Ménard et al., 2000).Since the SPREAD is tuned to match a subset of the elements of the state vector rather than the entire state vector, we do www.atmos-chem-phys.net/15/5835/2015/Atmos.Chem.Phys., 15, 5835-5850, 2015 Table 1.Results for the optimal runs (i.e., maximum wind extraction potential, WEP) for each experiment.The localization length (L) is provided along with WEP and global mean root mean square error (RMSE) for u, v, z, and q.NMI refers to normal mode initialization applied to the analysis fields.not expect χ 2 to be exactly 1, but it should be relatively close, at least in the time average.For these experiments, χ 2 (not shown) starts out slightly high, but levels out to a time mean (averaged from 2 to 10 days) of 0.99 for EnKF-uv and 0.97 for EnKF-ψχ .In Fig. 8, the analysis errors are projected onto the GW and RW modes.As expected from the imbalance calculations, the uninitialized EnKF-uv has much larger GW error due to larger imbalance in the increments.However, the EnKF-uv has slightly smaller RW errors.This is consistent with the initialized EnKF-uv having slightly larger WEP than the initialized EnKF-ψχ (Fig. 7a).This difference may be partly due to background error estimation biases caused by the ψχlocalization, as discussed in Kepert (2009).These biases will either overweigh or underweigh the background at different scales, resulting in a suboptimal solution.We could try to correct for this effect by altering the observation error covariance as in Kepert (2009), but this does not account for the scale dependence of the bias.The situation in our case is also complicated by the adaptive inflation, which uses different state variables in EnKF-uv and EnKF-ψχ.
Figure 9 (column 1) presents several global error diagnostics versus L for Experiment 1.We define the "optimal" runs for each experiment as those that maximize WEP.These are indicated by squares (triangles) for uninitialized (initialized) results in Fig. 9a and are also listed in Table 1.In Fig. 9a and b the uninitialized WEP and zonal wind errors (meridional wind errors are very similar and are not shown) exhibit strong dependence on L, with maximum WEP occurring at L = 1500 km for EnKF-uv and L = 2500 km for EnKF-ψχ.
The optimal EnKF-ψχ run results in larger WEP and smaller wind error compared EnKF-uv, which would appear to favor the choice of ψχ or uv.However, EnKF-ψχ is much more sensitive to variations in L, with WEP values actually becoming negative at small and large L.
Figure 9c shows a height error minimum at L = 2000 km for uninitialized EnKF-ψχ , while for uninitialized EnKFuv the height error increases monotonically with L. The increase of height error with L is driven by the increase of imbalance with L, as seen in Fig. 9d.This increase in imbalance with longer L when assimilating only ozone observations is a new finding, which is opposite to the tendency of localization to create imbalance when assimilating dynamical observations, as will be shown in Experiment 2 and discussed by Mitchell et al. (2002) and Kepert (2009Kepert ( , 2011)).It is likely that ozone observations cause increased imbalance with L due to spurious ensemble correlation between ozone and the dynamical variables at large distances, which are projected onto the gravity modes.Up to L = 4000 km, the EnKF-uv has larger imbalance than EnKF-ψχ , which is consistent with the single-observation simulations.The ozone errors (Fig. 9e) show a broad minimum, with the EnKF-uv providing slightly better results at all L values.
After application of NMI, for both EnKF-uv and EnKFψχ the wind and height errors are smaller and WEP is larger at all values of L. The ozone error does not change, since the NMI is applied only to the dynamical fields.The initialized results show WEP maximizing at ∼ 84 % for EnKF-uv and ∼ 81 % for EnKF-ψχ (see triangles in Fig. 9a).The length scales corresponding to these values increase to 2000 km for EnKF-uv and 3500 km for EnKF-ψχ , suggesting the correlations at larger lengths are more reliable.That EnKF-uv outperforms EnKF-ψχ when NMI is applied is consistent with Kepert (2009), who showed that EnKF-uv (with NMI) resulted in smaller height and wind errors than the EnKF-ψχ (with NMI) due to background error estimation biases caused by the ψχ localization.
Up to this point, we have examined globally averaged analysis errors.To determine regional impact, we also examine how the errors vary with latitude for Experiment 1, shown in Fig. 10 (column 1).The initial wind errors (black lines in Fig. 10a and b) are largest in the NH tropics and midlatitudes and near the North Pole (global maps of initial wind and height errors are provided in Fig. 3 of Allen et al., 2014).After assimilating ozone, the uninitialized wind errors are reduced at most latitudes.Small increases in uninitialized zonal wind error occur near 70 • S and 70 • N.That the tropical bias has been removed is important, since the tropical winds are not as well constrained in the stratosphere by radiance observations alone.The uninitialized height errors (Fig. 10c) are more uniform after ozone assimilation and show slight improvement in some regions.However, uninitialized height errors have also increased over large portions of the globe, particularly for EnKF-uv.This is due largely to the imbalance generated by the ozone observations.Results with NMI (dotted lines in Fig. 10a-c) show reduced height (and wind) errors at all latitudes compared to the original analyses, due to removal of spurious GW.Ozone errors (Fig. 10d) are also reduced at all latitudes in this ozone assimilation experiment, with slightly larger errors in the tropics.

Experiment 2: height only
We now examine the results of Experiment 2, when only height data are assimilated.For both EnKF-uv and EnKFψχ , Fig. 9 (column 2) shows that WEP initially increases and wind errors decrease with less localization (larger L), with optimal values occurring at L = 5000 km for EnKF-uv and L = 7000 km for EnKF-ψχ, followed by a slight degradation at larger lengths.The EnKF-ψχ generally performs better than EnKF-uv.The minimum wind errors are ∼ 1.3 m s −1 for EnKF-uv and 1.0 m s −1 for EnKF-ψχ , which are reasonable values for a well-constrained stratospheric analysis.For example, Hertzog et al. (2004) compared NH stratospheric analyses with observations from long-duration balloon flights and calculated error SDs of the zonal wind components of ∼ 1.3 m s −1 for ECMWF and ∼ 1.9 m s −1 for NCEP, when the observations were low-pass filtered to remove the variance due to inertia-gravity waves.The biggest difference between experiments 1 and 2 is that assimilation of height observations results in much less imbalance (note different vertical scales in Fig. 9d and i).The imbalance, like the wind and height errors, generally decreases with less localization, which is opposite to what occurred in Experiment 1.However, this is consistent with previous studies that have examined balance in the EnKF in the context of assimilation of dynamical variables (e.g., Mitchell et al., 2002).This result provides a caution that while reducing the localization may reduce imbalance for some observations, it may increase imbalance when assimilating ozone.Applying NMI to the analyses results in almost no change to the WEP and wind errors, but does improve the height errors, particularly for EnKF-uv.
The errors as a function of latitude for Experiment 2 are shown in Fig. 10 (column 2).The wind errors (Fig. 10e,  f) are largest in the tropics and decrease towards the poles.This is expected, since the height is more strongly correlated with wind in the extratropics due to geostrophic balance.In the tropics this balance breaks down, and it is more difficult for the EnKF to correct the winds with height observations alone.The analyzed height errors (Fig. 10g) are markedly reduced from the initial errors, with slightly larger values in the extratropics.Experiment 2 also improves the ozone, but only by a small amount (Fig. 10h).The small ozone improvement is in the extratropics, likely due to more accurate winds that drive ozone advection.

Experiment 3: ozone and height
The final experiment examines the value of adding ozone assimilation to the analyses produced by the height only assimilation.The results as a function of L are shown in Fig. 9 (column 3).This experiment results in the smallest errors and highest WEP values, confirming that ozone and height observations provide complimentary information to the DA system.Large WEP values occur for a broad range of L, indicating that the results are not very sensitive to the choice of localization length.The lowest uninitialized wind errors occur at L = 3500 km for both EnKF-uv and EnKF-ψχ , and the maximum uninitialized WEP (∼ 87 %) is larger than when either ozone or height are assimilated separately.The application of NMI slightly increases the optimal WEP to ∼ 90 % for both systems.
The uninitialized height error and imbalance for Experiment 3 (Fig. 9m, n) show broad minima, which reflects the combined tendencies of the height observations to increase imbalance at small L values and the ozone observations to increase imbalance at large L values.The imbalance remains relatively low in these experiments (< 20 m for L < 5000 km), with the EnKF-uv showing somewhat higher values than EnKF-ψχ .The error in initialized height (dashed lines of Fig. 9m) does show a significant decrease, suggesting that there is some GW noise in Experiment 3, but it is much less than when ozone is assimilated alone.It appears that combining height observations with ozone observations reduces the GW that would otherwise be generated by the ozone observations alone.
The errors as a function of latitude for the optimal results from Experiment 3 are presented in Fig. 10 (column 3).Wind errors are quite small (< 0.5 m s −1 ) across the globe, with the primary benefit occurring in the tropics, where wind errors are much less than either Experiment 1 or 2. The ozone in Experiment 3 is also better than in Experiment 1. Having a  and q [ppmv] (rows 1, 2, 3, and 4, respectively) for the optimal runs (as shown in Table 1 and in the highlighted squares of Fig. 9) of experiments 1, 2, and 3 (columns 1, 2, and 3, respectively).Black lines show initial errors and red (blue) lines show EnKF-uv (EnKF-ψχ) errors.Solid (dotted) lines indicate uninitialized (initialized) results (there are no dotted lines in row 4 because the ozone error does not change, since the NMI is applied only to the dynamical fields).
better background ozone field (due to better winds) allows for more efficient use of the ozone observations.This likely provides a positive feedback in the system that enhances the ozone impact.The addition of ozone also tends to flatten height errors with latitude.Application of NMI does not impact the wind errors very much, but does reduce the height errors as seen in Fig. 10k.
Additional experiments (not shown) were performed with a much smaller height error SD of 10 m.Decreasing the height error for height only assimilation increased the maximum WEP to 84 % for EnKF-uv and 89 % for EnKF-ψχ .The impact of adding ozone assimilation in these experiments was also positive, increasing the maximum WEP to 92 % for EnKF-uv and 91 % for EnKF-ψχ , with most of the wind improvements occurring in the tropics.These results suggest that even in a very well-tuned system, high-quality ozone observations can, in principle, improve the winds.The EnKF DA is able to employ cross-correlations between state variables in the ensemble background states to couple tracer and dynamical variables.This study examined several aspects of extraction of wind information from EnKF ozone assimilation using a SWM coupled with ozone advection.Three sets of experiments were performed that assimilated ozone, geopotential height, or both.Modest improvements to the winds were observed when either ozone or height were assimilated separately.Final WEP values of 46 % (58 %) were obtained for ozone and 60 % (69 %) for height with the EnKF-uv (EnKF-ψχ ) system.When NMI was applied to the ozone experiment, WEP jumped to 84 % (81 %), showing that gravity wave noise was generating significant error.The NMI applied to the height experiment resulted in WEP increases of less than 1 %, suggesting very small imbalance.When assimilating both ozone and height, WEP rose to ∼ 87 % for both systems.Imbalance was also much less than when ozone was assimilated alone.The addition of height observations appears to reduce the gravity wave noise in the EnKF DA, thereby reducing the need for initialization.This is important, since over-filtering could be a problem if NMI is applied to the upper stratosphere/mesosphere (Sankey et al., 2007) and the tropics (Nezlin et al., 2009), where unbalanced modes play an important role in the real atmosphere (see also Koshyk et al., 1999).Applying NMI to the combined experiment resulted in a modest increase in WEP to ∼ 90 %.The greatest impact of ozone assimilation on the winds was found to occur in the tropics, which are less well constrained by height assimilation due to lack of geostrophy.This study also compared results from EnKF systems that used different flow variables.While the EnKF-uv system with ozone observations generated greater imbalance, this system was also able to more accurately determine the wind structure of the rotational wave modes.This may be due to biases in the specification of the background error covariance in the ENKF-ψχ, as discussed by Kepert (2009).As a result, when NMI was applied, the EnKF-uv performed slightly better than the EnKF-ψχ .For height assimilation, the EnKFψχ performed better, due to less imbalance, while the combined assimilation of ozone and height produced similar results in the two systems.
In each experiment the localization length was tuned to maximize the wind extraction.Previous studies have shown that tighter localization increases imbalance, which may be detrimental.We showed that while this was the case for height observations, for ozone observations the imbalance actually increased with localization length.The cause is uncertain, but may be due to spurious long-range correlations between ozone and the dynamical fields, which project onto the gravity modes of the SWM.
While under the ideal conditions used in this study WEP values of up to ∼ 90 % were achieved (wind component er-rors ∼ 0.3 m s −1 ), there are many challenges to demonstrating that the ozone-wind coupling in an operational DA system can be beneficial.There are observation system challenges such as frequency, latency, precision, and bias.There are also modeling challenges such as accurate ozone transport, chemistry, and radiation.The results here were obtained with a single-layer model, relatively low resolution (T42), and a rather simple wave-forcing scenario.Given these caveats, this study demonstrated ozone-wind interaction in the EnKF and the potential for ozone assimilation to benefit the wind analysis, particularly in the tropics.
Whether 4D-Var or EnKF is better for ozone-wind extraction is still an open question.In our previous work we showed that wind extraction is feasible when assimilating globally gridded hourly tracer data (ozone, nitrous oxide, or water vapor) within 4D-Var.Follow-up experiments (not presented here) indicate that ozone-wind extraction is also possible in 4D-Var assimilation of the ozone and height data used here.In future work, we plan to directly compare 4D-Var, EnKF, and hybrid methods for tracer-wind extraction.

Figure 3 .
Figure 3.The normal mode frequencies [day −1 ] for the spectral SWM at triangular truncation T42 as a function of zonal wave number (m) for the first six values of the meridional wave number l = n− |m|, where n is the total wave number.Positive (negative) values of m indicate eastward (westward) motion.Modes are separated into (a) eastward and westward gravity wave (GW) modes, and (b) rotational wave (RW) modes.There are no rotational modes for positive m (see Sect. 2.3).The cutoff frequency (1.0 day −1 ) for the NMI is indicated by the dotted line.Note that the frequency scales are different for the two plots for easier viewing.

Figure 6 .
Figure 6.Sampling patterns for 24 h of (a) polar-orbiting ozone data and (b) pseudo-random height data.

Figure 7 .
Figure 7. Diagnostics from optimal runs of Experiment 1: ozone only (L = 1500 km for EnKF-uv and L = 2500 km for EnKF-ψχ ).(a) WEP [%], (b), (c), and (d) RMSEs for u [m s −1 ], v [m s −1 ], and z [m], respectively, (e) imbalance [m], and (f) RMSE for q [ppmv].EnKF-uv is red and EnKF-ψχ is blue.Solid (dotted) lines indicate uninitialized (initialized) results (there are no dotted lines in (f) because the ozone error does not change, since the NMI is applied only to the dynamical fields).In (f) the ozone observation error standard deviation is indicated by the horizontal dotted line.Blue circles at day 0 indicate the initial values.

Figure 10 .
Figure 10.RMSEs as a function of latitude for u[m s −1 ], v [m s −1 ], z [m],and q [ppmv] (rows 1, 2, 3, and 4, respectively) for the optimal runs (as shown in Table1and in the highlighted squares of Fig.9) of experiments 1, 2, and 3 (columns 1, 2, and 3, respectively).Black lines show initial errors and red (blue) lines show EnKF-uv (EnKF-ψχ) errors.Solid (dotted) lines indicate uninitialized (initialized) results (there are no dotted lines in row 4 because the ozone error does not change, since the NMI is applied only to the dynamical fields).