Wind extraction potential from 4 D-Var assimilation of stratospheric O 3 , N 2 O , and H 2 O using a global shallow water model

The wind extraction due to assimilation of stratospheric trace gas (tracer) data is examined using a 4DVar (four-dimensional variational) data assimilation system based on the shallow water equations coupled to the tracer continuity equation. The procedure is outlined as follows. First, a nature run is created, simulating middle stratospheric winter conditions. Second, ozone (O 3), nitrous oxide (N2O), and water vapor (H2O) (treated in this study as passive tracers) are initialized using Aura Microwave Limb Sounder (MLS) mixing ratios at 850 K potential temperature and are advected by the nature run winds. Third, the initial dynamical conditions are perturbed by using a 6 h offset. Fourth, simulated hourly tracer observations on the full model grid are assimilated with a 4D-Var system in which tracer and winds are coupled via the adjoint of the tracer continuity equation. Multiple assimilation experiments are performed by varying the amount of random observation error added to the simulated measurements. Finally, the wind extraction potential (WEP) is calculated as the reduction of the vector wind root mean square error (RMSE) relative to the maximum possible reduction. For a single 6 h assimilation cycle with the smallest observation error, WEP values are ∼ 60 % for all three tracers, while 10-day multi-cycle simulations result in WEP of∼ 90 %, wind errors of∼ 0.3 m s−1, and height errors of∼ 13 m. There is therefore sufficient information in the tracer fields to almost completely constrain the dynamics, even without direct assimilation of dynamical information. When realistic observation error is added (based on MLS precisions at 10 hPa), the WEP after 10 days is 90 % for O3, 87 % for N2O, and 72 % for H2O. O3 and N2O provide more wind information than H 2O due to stronger background gradients relative to the MLS precisions. The RMSE for wind reach a minimum level of ∼ 0.3–0.9 m s −1 for the MLS precisions, suggesting a limit to which realistic tracers could constrain the winds, given complete global cover age. With higher observation noise levels, the WEP values decrease, but the impact on the winds is still positive up to noise levels of 100 % (relative to the global mean value) when compared to the case of no data assimilation.


Introduction
This study examines the problem of extraction of global wind information from the assimilation of stratospheric trace gas (tracer) observations. The assimilation of stratospheric ozone in 4D-Var (four-dimensional variational) algorithms and the resultant dynamical coupling have previously caused problems in operational numerical weather prediction (NWP) (Han and McNally, 2010;Dragani and McNally, 2013). Among the tracer assimilation challenges are paucity of observations, insufficient data quality, insufficient modeling of tracers in the forecast model, as well as phenomenological limitations due to tracer variations with respect to the wind vector. The problem of tracer-wind extraction has been examined using a series of increasingly sophisticated models ranging from 1-D (Daley, 1995;Allen et al., 2013), 2-D (Daley, 1996;Riishøjgaard, 1996), and 3-D (Peuch et al., 2000;Semane et al., 2009;Allen et al., 2013). A review of the results from these studies is provided in Allen et al. (2013). The basic conclusion is that with complete data coverage, sufficient precision, a high-quality background tracer field, and a 4D-Var or extended Kalman filter algorithm, wind extraction is possible across the globe and throughout the stratosphere and upper troposphere. However, the ability to extract wind is greatly hindered by limitations in any of these parameters, making the current value of operational stratospheric tracer-wind coupling in NWP systems questionable.
Given the observational and modeling challenges, one may ask whether the continued pursuit of stratospheric tracerwind extraction is worthwhile. One main motivation is the current deficiency of adequate global wind profile information (World Meteorological Organization, 2000). Upper air wind observations from pilot reports, radiosondes, and cloud and water vapor feature-tracking leave large gaps, particularly in the tropics, Southern Ocean, and in most of the stratosphere and mesosphere. In the extratropics, traditional mass-wind balance provides some constraint on the wind from temperature information based on radiance observations, but this balance fails in the tropics and for sub-synoptic scales (less than ∼ 1000 km) in the extratropics. While spaceborne Doppler wind lidar (DWL) has the potential to provide wind profile measurements in the troposphere and lower stratosphere, it is extremely expensive (National Research Council, 2007), and the only currently planned DWL mission, Atmospheric Dynamics Mission (ADM-Aeolus) (Stoffelen et al., 2005), will be limited to a single line-of-sight wind component, altitudes below ∼ 26 km, and simple alongtrack sampling. Another motivation is that 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 show that the dominant factor is the adjustment of the wind field in order to match observed humidity features (the so-called "tracer advection effect").
While the tropospheric benefit of tracer-wind interaction has been established, assimilation of stratospheric tracer observations may also provide stratospheric wind information (Semane et al., 2009;Milewski and Bourqui, 2011;Allen et al., 2013). The studies to date on stratospheric trace gas assimilation have focused on ozone (O 3 ) due to its relatively long lifetime, the roles it plays in the radiative heating/cooling in NWP systems and infrared radiance assimilation, and the availability of high quality satellite observations (e.g., from the Suomi National Polar-orbiting Partnership -NPP -Ozone Mapping and Profiler Suite -OMPS; Flynn et al., 2009). There are other long-lived stratospheric tracers, however, that could be considered as potentially useful for wind extraction. These include nitrous oxide (N 2 O), water vapor (H 2 O), methane (CH 4 ), and chlorofluorocarbons. Modeling each of these species in an NWP system poses many challenges, but in principle these tracers contain dynamical information that could be "harvested" to benefit the stratospheric winds. This paper examines the relative impact on the winds due to assimilation of three of these tracers: ozone, nitrous oxide, and water vapor.
The experiments are performed using a 4D-Var data assimilation system based on the shallow water model (SWM) equations. Application of the SWM equations to the variational assimilation problem has proved useful for algorithm development and basic understanding of the problem (e.g., Courtier and Talagrand, 1990). The SWM system provides sufficient complexity that realistic flow situations can be simulated, while maintaining greater control of the system parameters than in the full NWP system. This study builds on the 2-D tracer assimilation studies by Daley (1996) and Riishøjgaard (1996). In the stratosphere, which is the focus of this paper, the tracer-wind extraction should be to first order a 2-D quasi-isentropic problem, since the horizontal transport timescales are smaller than the vertical timescales (Plumb, 2002). Near the tropopause, strong vertical ozone gradients and vertical motion make the full 3-D modeling necessary (Riishøjgaard, 1996). The experiments in this study use an observation grid that matches the model grid and provides complete global coverage in order to maximize wind extraction. The goal is to determine the relative wind extraction potential due to the geophysical variability of each of the tracers, rather than to determine the impact of the observing grid. This idealized approach provides theoretical estimates of the amount of dynamical information that can be extracted from each tracer. Future work is planned to examine the impact of various observation systems as well as simultaneous assimilation of tracer and dynamical fields.
Section 2 describes the SWM data assimilation system. Section 3 describes the experimental design, and Sect. 4 presents the results of both single-cycle and multiple-cycle experiments, focusing on the relative impact of the different tracers. Sections 5 and 6 provide a discussion and conclusions, respectively.

Forecast model
The forecast model is a spectral model based on the vorticity divergence form of the SWM equations as formulated in Sect. 2a of Ritchie (1988), with the inclusion of fourth-order diffusion applied to the vorticity, divergence, and geopotential following the semi-implicit formulation. A tracer advection equation is added, which solves for the mixing ratio of a passive tracer as a function of time with the same fourthorder diffusion operator. The system is run at triangular truncation T42 with a time step of 1200 s for both the forward model and for the tangent linear model (TLM), and adjoint integrations discussed in Sect. 2.2. The diffusion coefficient is set to 5.0 × 10 15 m 4 s −1 , as in Jakob-Chien et al. (1995). The model fields are saved on the Gaussian grid with 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). The global mean geopotential height (i.e., the free upper surface of the fluid) is specified to be 10 km. The algorithm is very similar to the model described by Jakob-Chien et al. (1995). The forecast model was validated using tests developed by Williamson et al. (1992) and Galewsky et al. (2004), with excellent results.

Data assimilation system
We employ a 4D-Var data assimilation system (DAS) using the accelerated representer (AR) approach (Xu et al., 2005;Rosmond and Xu, 2006), which is an observation space algorithm based on the solution to a standard 4D-Var cost function defined as follows (note: the "perfect model" assumption is used, although not required for the AR approach).
Here, J b 0 is the cost function for deviations of the initial analyzed state vector x 0 from the initial background state x b 0 . In the DAS algorithm, the horizontal flow is converted from the forecast model variables (vorticity and divergence) to zonal and meridional wind (u and v). Although either combination of variables could be used, the background errors are more straightforward to model with the observable variables u and v. P b 0 is the initial background error covariance matrix. J r is the cost function for the observations, where y n is the vector of observations at time t n , H (x n ) is the nonlinear observation operator that maps the model state vector onto the observations, and R is the observation error covariance matrix. The cost function is minimized over a 6 h analysis window comprising N = 18 time steps. The representer method is used to obtain the solution, incorporating the linearized forecast model M and the linearized observation operator H: and d = y − H (x b ) is the innovation vector. The background state is calculated by integrating the nonlinear forecast model M, subject to the initial conditions x b 0 .
The analysis is obtained by first solving [H P b H T + R] z = d using a standard conjugate gradient descent algorithm to obtain z = [H P b H T + R] −1 d and then post-multiplying Table 1. Schematic showing the matrix of error covariances that are modeled in P b 0 . σ is the background error standard deviation, and ρ is the two-point spatial correlation between variables. The rows and columns are ordered by state variable as follows: row/column 1: u (zonal wind), row/column 2: v (meridional wind), row/column 3: (geopotential), and row/column 4: q (tracer). Each table entry represents the block covariance matrix for the two specified variables. Zeros indicate no correlation for these two variables. Figure B2 of Daley and Barker (2001, available online) shows an example of the spatial structure of the correlations.
(2), both the TLM, M, and its transpose (or adjoint), M T , are integrated over the analysis window during each iteration of the descent algorithm in order to "sweep through" all observations and to calculate the time-evolving background error covariance P b . In 4D-Var, the coupling of tracer and dynamical fields (wind and geopotential) occurs implicitly via the adjoint of the TLM, which propagates sensitivities of the cost function with respect to the tracer observations backwards in time to obtain sensitivities of the cost function to the initial model state (see Allen et al., 2013, for further discussion). The solution is iterated until the gradient of the cost function drops by a factor of 20 (ratio to initial gradient of 0.05) (see Daley and Barker, 2001 for details). If the solution does not reach this threshold, the solver is terminated after 100 iterations. The TLM is constructed from the full nonlinear model by expanding each forecast model state variable into time-evolving background and perturbation terms and removing terms that are quadratic in perturbations. The quality of the TLM, as discussed in Sect. 3, is quite good over the 6 h assimilation window for the initial conditions used in this paper. The adjoint model was constructed line-by-line from the TLM following the adjoint algorithms of Giering and Kaminski (1998) and guidance from Rosmond (1997). Each piece of TLM and adjoint code was validated using a standard adjoint operator test to ensure correct coding.
The initial background error covariance, P b 0 , includes multivariate correlations based on approximate geostrophic balance. Table 1 illustrates the covariance elements calculated for any pair of grid points. The nine correlations between the dynamical variables (zonal wind, meridional wind, and geopotential) are calculated using the equations in Appendix B of Daley and Barker (2001), which are spherical extensions of the tangential plane equations in polar coordinates derived from Chapter 5 of Daley (1991). The wind-geopotential correlations are based on approximate geostrophic balance, with a coupling between the geopotential and the stream function that varies with latitude as illustrated in Fig. 1b. The wind-wind correlations are derived from the stream function correlations using the Helmholtz equations. The spatial correlation of the stream function is isotropic and homogeneous and is modeled as the single-order auto regressive (SOAR) function with a prescribed horizontal correlation length. Various values of the horizontal correlation length were tested to maximize tracer-wind extraction. We decided to use a length that varied smoothly with latitude from 1500 km in the tropics to 1000 km in the extratropics (see Fig. 1a). This length scale works well for the simulations done here at T42, but may need to be altered for simulations at higher resolution. The ratio of the divergent kinetic energy to total horizontal kinetic is specified as 0.05 (see Chapter 5 of Daley, 1991 for details). There is no coupling between the tracer and dynamical variables through P b 0 , as illustrated by the zeros in Table 1. There are also no spatial correlations in the P b 0 for the tracer field, thus allowing tracer increments of very fine scale. The entire P b 0 matrix is too large to store in memory, so it is calculated during each iteration of the solver, with several terms in P b 0 precomputed for efficiency.
In order to optimize the wind extraction, the background error standard deviations (σ b ) for zonal wind (σ u ), meridional wind (σ v ), and tracer (σ q ) are updated for each assimilation cycle. They are calculated as the globally averaged standard deviation of the difference between the previous analysis and the truth. This allows the background errors to decrease as the analysis improves, a self-tuning approach that works well in these idealized experiments. The initial σ u is 3.2 m s −1 (see Fig. 1c), which represents the globally averaged 6 h wind differences from the nature run, while the initial σ q is zero. The background error standard deviation for geopotential (σ ) is calculated from the wind values using the cross-correlation functions discussed above. To avoid small background errors where geostrophic balance breaks down, σ is fixed in the tropics at the value calculated at 30 • latitude. Figure 1d shows the initial σ for these settings. While the background error standard deviations vary over the 10-day experiments, the correlation length scales do not change, and hence all the correlation functions (ρ uu , ρ vv , . . . ) remain fixed with time.
Because the observations are defined on the Gaussian grid, no interpolation is necessary. The observation operator H is trivial, and H = H . R was defined as diagonal for these experiments, with no spatial or temporal correlations between observations. The observation error standard deviation σ ob is a constant value in terms of the tracer mixing ratio, and will be specified below.

Experimental design
The wind extraction experiments follow the design of Allen et al. (2013) in which a nature run is used as the truth, both for the tracer and dynamical fields (winds and geopotential). This run is then used to create "perfect" observations on the full T42 Gaussian grid. Gaussian random errors are added to the perfect observations, as specified for each experiment. An initial, imperfect, 6 h background forecast (or "background trajectory") is then created by adding perturbations to the dynamical fields of the nature run. A 6 h data assimilation cycle is performed to create a 6 h "analysis trajectory". Comparing the background and the analysis with the nature run allows determination of the impact of tracer assimilation on the dynamical fields. Below we describe each component of these wind extraction experiments in more detail. We note here that since the same forecast model is used for the nature run and for the assimilation, the results are likely over-optimistic. However, the intent of this study is not to determine the impact of any realistic observing system, but rather to examine the theoretical limitations of the wind information contained within the tracer field.
The nature run was designed to simulate Northern Hemisphere (NH) winter conditions in the middle stratosphere. This was based on previous simulations by Juckes (1989( ), Norton (1994, and Thuburn and Lagneau (1999), in which an SWM system is forced by the bottom topography being raised and lowered in order to simulate planetaryscale waves. The model was initialized with 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). The jet strength is set by u max = 60 m s −1 . The initial geopotential height h = /g, where is geopotential and g is the Earth's gravitational acceleration, is specified using the gradient wind balance and given a mean height of 10 km. The topographic forcing is produced following the approach of Norton (1994). A mountain of height 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 starts at day 20 of this nature run, there is therefore no surface topography during the assimilation.
In the nature run the tracer is initialized using the Aura Microwave Limb Sounder (MLS) trace gas 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. Simulations were performed with O 3 , N 2 O, and H 2 O. The tracer in each case is treated as passive (i.e., no chemical source/sink) and there is no radiative interaction between tracer and dynamics. N 2 O and H 2 O have relatively long photochemical lifetimes in the middle stratosphere, as does O 3 in the winter hemisphere. However, the photochemical lifetime of O 3 becomes shorter (about 10 days) in the tropics and summer hemisphere at 10 hPa (e.g., McCormack et al., 2006, Fig. 2), and therefore photochemical affects will become more important. The impacts of chemistry and radiation on the wind extraction process are reserved for future work. Figure 2 plots the potential vorticity (PV = ζ a / h, where ζ a is the absolute vorticity) in 5-day intervals, starting on the 15th day during the nature run along with the three tracers. The wave forcing produces an "Aleutian high" (centered at 30 • longitude) type wave response in the NH with the polar vortex (centered near 180 • longitude) pushed off the pole. Low O 3 and N 2 O (high H 2 O and PV) are visible in the polar vortex and high O 3 and N 2 O (low H 2 O and PV) are entrained into the "Aleutian high". The sharpening of the PV gradients with time along the vortex edge is seen, as expected for these types of wave-breaking simulations. The three tracer fields have very similar contour shapes, as expected for long-lived tracers in a 2-D flow (Rhines and Young, 1983). However, the dynamic range (ratio of maximum to minimum value) of the tracers is considerably different, resulting in different horizontal gradients. The dynamic ranges at initialization (day 0) of the nature run are as follows: O 3 [11.0/4. To start the assimilation experiments, the background forecast is initialized from the nature run using the true tracer distribution but with dynamical fields sampled at a 6 h offset. This method allows for perturbed initial dynamical fields that are approximately in balance. Although the tracer is perfect at the first time step, tracer errors develop over the sixhour background forecast by advection with imperfect dynamical fields. The first assimilation cycle represents the best case scenario for wind extraction because all the tracer error, and hence the deviation relative to simulated observations, is caused by wind error. On further cycling of the assimilation, the background tracer field will retain error from the previous cycle in addition to new error from imperfect advection. The assimilation is initialized 20 days into the nature run (T = 20 d, 0 h), and the initial dynamical fields are then taken from day 20 + 6 h, while the tracer field is taken from day 20 + 0 h. The background forecasts always use the full forecast model with the same settings as the nature run. In the multi-cycle simulations, the background for each successive cycle is produced from an integration starting with the previous 6 h analysis.
Observations are simulated at each hour and at each Gaussian grid point by sampling the nature run and adding varying amounts of random error based on percentage of the global mean mixing ratio (mean values are 8.375 ppmv (parts per million by volume) for O 3 , 0.137 ppmv for N 2 O, and 5.219 ppmv for H 2 O). For each experiment, the σ ob used in the 4D-Var solution was specified consistently with the applied random error, so we use the symbol σ ob to represent both terms. Ten-day experiments were performed for σ ob of 0.1-1 % at 0.1 % increments, 1-10 % at 1 % increments, and 10-100 % at 10 % increments. When σ ob becomes sufficiently small, the 4D-Var cost function minimization algorithm fails to converge. For each tracer, the lowest σ ob error for which the algorithm converges in 100 iterations is labeled as the base error standard deviation (BEST) experiment. The BEST experiments correspond to σ ob of 0.5 % for O 3 , 2 % for N 2 O, and 0.3 % for H 2 O. As will be shown (Fig. 10), as the observation errors decrease, wind extraction asymptotically approaches a maximum, such that the exact choice of BEST observation error is not critical. In addition, we performed an experiment in which we specified σ ob consistent with MLS retrieval error at 10 hPa (∼ 30 km), with values of 0.1 ppmv for O 3 , 0.01 ppmv for N 2 O, and 0.3 ppmv for H 2 O (see MLS Version 3.3 level 2 data quality and description document, Livesey et al., 2011). These are referred to as MLS experiments. The MLS values are used here to represent typical state-of-the-art measurements, even though the simulated observation density and frequency are better than any foreseen measurement system. Note that σ ob for all experiments is in mixing ratio units, and the relative error will vary depending on the tracer distribution.
As a global consistency check of the solution, we calculated , which should equal the number of observations (Ménard et al., 2000). These two numbers agree to within 1.5 % (based on root mean square error (RMSE) differences over 10 days) for all experiments shown.

Results
We first examine the results of the first analysis time step, or single-cycle assimilation. This allows us to compare wind extraction potential for each tracer under the same initial conditions (identical background wind error and perfect initial background tracer fields). We then examine results for 10 days of tracer assimilation to determine how well tracer assimilation alone is able to constrain the model dynamics. In a typical NWP system, tracer data is not assimilated alone, but in combination with a full suite of meteorological observations. However, the goal of this study is to examine the limits of wind extraction solely from tracer measurements.
First, we examine the BEST O 3 assimilation. Figure 3 shows the zonal average background and analysis errors for the first analysis cycle. In order to verify the analysis over the entire time window, the analysis error is plotted for 0, 3, and 6 h. Global maps of the analysis increments (analysis -background) and ideal analysis increments (truth -background) at 0 h are also shown for each of the dynamical variables. For all three variables, the analysis errors are less than the background errors over most of the globe. The most dramatic error reduction occurs for the meridional wind, where errors drop from ∼ 2-4 m s −1 down to ∼ 1.0-1.5 m s −1 . For geopotential, the reduction is greatest in the extratropics of both hemispheres. Winds are affected by tracer assimilation directly via the adjoint of the tracer continuity equation. The influence on the geopotential is secondary via advection of the geopotential and via the wind-geopotential correlations in P b 0 . Since the latter is absent in the tropics (see Fig. 1b), it is expected that the tropical height field will be less influenced by tracer assimilation. The wind analysis error is similar throughout the analysis time window. The height analysis error, however, shows more sensitivity with time, with smaller RMSE in the extratropical SH and the subtropical NH toward the end of the analysis window. In the NH extratropics, the height errors actually increase with time over the analysis window. The ability to extract wind is partially dependent on the accuracy of the TLM integration compared to the full nonlinear model. The TLM error is also plotted in Fig. 3. Non-zero errors indicate the effects of neglecting the nonlinear terms in the SWM equations. The TLM errors are everywhere smaller than the analysis errors. The zonal wind TLM error shows a peak at 40 • N, suggesting that nonlinear terms are playing a role there. The TLM errors give an estimate of the limit of 4D-Var wind extraction for this first cycle. The zonal average TLM wind errors are less than ∼ 0.5 m s −1 across the globe and height errors are less than ∼ 5 m, and both are clearly a small percentage of the analysis errors. Comparing the increment maps in Fig. 3, we see that most of the large-scale features in the analyzed wind and height increments match those in the ideal increments. The increments are generally larger in the NH, which is expected since the dynamical fields of the nature run were forced by planetary waves in the NH.
We next compare BEST single-cycle results for all three tracers. Figure 4 (top row) shows the background and analysis errors for O 3 , N 2 O, and H 2 O. The error reductions in zonal and meridional wind for all three tracers are similar. There are more noticeable differences in the height errors, particularly in the NH extratropics. The single-cycle results for experiments using stratospheric MLS precision estimates for σ ob are also shown in Fig. 4. Unlike the BEST results, the wind and height error reduction varies by tracer, with O 3 having the strongest impact, followed closely by N 2 O. H 2 O assimilation in this experiment appears to have much less of an impact on the dynamical fields.
We further quantify these results by introducing the wind extraction potential (WEP), which we define as follows. First, we calculate the analyzed RMSE for zonal and meridional wind errors as a function of latitude and time, u RMSE (φ, t) and v RMSE (φ, t), and from these we obtain the vector wind error The percentage difference in vector wind error relative to the initial error is then calculated, × 100 %, and WEP is defined as the area-weighted global mean of this quantity. A WEP value of 100 % indicates the analysis equals the truth (V RMSE (φ, t) = 0). Table 2 lists the WEP values for select experiments. For BEST, the WEP after one cycle is similar for all three tracers (∼ 57-59 %), indicating that tracer assimilation alone is able to correct for almost 60 % of the wind error in only one cycle. For the MLS experiment, the WEP varies by tracer with O 3 having the strongest impact, and nearly the BEST impact, (57 %), followed closely by N 2 O (45 %). H 2 O assimilation results in much smaller error reductions (16 % WEP) than the other tracers. These results can be qualitatively explained by dividing the global tracer range (maximum -minimum mixing ratio) by the precision of each tracer (the range-to-precision ratio). This ratio (at 0 h of the first assimilation cycle) is ∼ 73 % for O 3 , ∼ 23 % for N 2 O, and ∼ 6 % for H 2 O. Although the relative precision value for N 2 O is larger than that of H 2 O, the dynamic range of N 2 O values is greater, causing it to give a larger WEP value. This is a simple rule-of-thumb relationship. More precisely, the wind extraction is influenced by the magnitude and orientation of the local tracer gradients relative to the wind errors (see Allen et al., 2013, for further discussion).
These simulations approximate only one stratospheric level (850 K or ∼ 10 hPa) and one set of dynamical conditions. The results for tracer distributions and dynamics for other levels and conditions will result in varying WEP. For example, we ran experiments using a nature run initialized with MLS data at 460 K (∼ 70 hPa) and used the MLS precision values of that level for σ ob , but used the same dynamical fields. The resulting WEP values after one cycle (45 % for O 3 , 20 % for N 2 O, and 3.2 % for H 2 O) were smaller than at 850 K. This can be qualitatively explained by smaller rangeto-precision ratios at 460 K (∼ 37 for O 3 , ∼ 12 for N 2 O, and ∼ 4 for H 2 O) relative to those at 850 K. While the precision value of O 3 is smaller at 460 K (0.04 ppmv), the range is much smaller due to lower mixing ratios. For N 2 O, the precision value is larger at 460 K (0.02 ppmv), while the range is slightly smaller. For H 2 O, the precision is the same as at 850 K (0.3 ppmv), but the range is smaller. These results are consistent with Allen et al. (2013), who showed greater wind extraction in the middle and upper stratosphere than in the lower stratosphere for a 4D-Var ozone assimilation experiment using a full NWP system. We next examine the results of 10 days of assimilation, corresponding to days 20 to 30 of the nature run. Although the initial background tracer field is perfect for the first cycle, background tracer errors develop over subsequent cycles due to imperfect wind fields. For these runs, we continue to use observation coverage on the full Gaussian grid at hourly intervals. Figure 5 shows the time variation of the globally averaged RMSE analysis for the BEST experiments. Like the single-cycle results, all three tracers show similar error reduction. By day 6, the zonal and meridional wind errors reach a minimum of ∼ 0.3 m s −1 , and the height errors drop to ∼ 12-13 m. The tracer errors, scaled by σ ob , are also shown in Fig. 5. The tracer errors increase over the first cycle from initially perfect conditions, and then decrease steadily, with tracer errors of < 15 % after 10 days. Figure 6 shows the zonal average, day 10 analysis errors for the BEST experiment as a function of latitude. Wind errors are nearly constant in latitude, with zonal wind error slightly larger in the NH. Height errors are also nearly constant, except near the poles. Tracer errors are generally larger in the tropics. Some latitude dependence may be expected due to the sampling of tracer observations on the Gaussian grid. Because there are no spatial correlations in the simulated observations or background covariance, the denser coverage near the poles should, in principle, preferentially reduce tracer errors near the poles.
In Fig. 7, the time variation of the analysis errors for the MLS experiments is shown. The winds are less constrained than in the BEST experiments with final values of  Table 2. Specifications for select assimilation experiments used in this paper. Wind extraction potential (WEP) is percentage of potential global vector wind reduction (see text for details). For the WEP after 10 days (last column), the value shown is the average value over the past four cycles. Note that ppmv refers to parts per million by volume. The WEP results over 10 days of assimilation are shown in Fig. 9 for BEST and MLS experiments, along with the 100% (relative to global mean) error experiment. For the  and "no data" (black). The "no data" indicates the WEP diagnostic applied to a 10-day forecast with no data assimilated. Circles indicate the initial WEP values (0 % in each case).
BEST experiment, the WEP value after 10 days is ∼ 90 % for all three tracers (see Table 2 for summary of WEP values). It is clear that the high precision and global hourly observation sampling in these experiments is sufficient to significantly reduce the initial wind errors. With errors based on MLS precision, the relative impact on winds from the three tracers is consistent with the single-cycle results, with day 10 WEP values of 90 % for O 3 , 87 % for N 2 O, and 72 % for H 2 O. For the 100 % error experiment, day 10 WEP is 12 % for O 3 , 36 % for N 2 O, and −19 % for H 2 O. Note that the negative WEP for H 2 O does not necessarily indicate a negative impact of assimilation. A calculation of day 10 wind errors resulting from a "no data" experiment (i.e., a 10-day forecast) results in a day 10 WEP of −30 %, and the time series of WEP for the "no data" experiment is always below the WEP with assimilation (see black line on Fig. 9). So assimilation of the tracers results in positive impact on the winds even for the case of 100 % relative error.
The WEP values obtained for the 10-day cycling experiments for BEST and MLS are very similar for O 3 and N 2 O, which suggests that the MLS precision for these tracers provides nearly the maximum possible wind extraction potential. To further examine this, we performed assimilation experiments using σ ob ranging from the BEST values up to 100 % of the global mean. The resulting WEP values are plotted in Fig. 10 as a function of σ ob normalized by the tracer global mean. For each tracer, the WEP asymptotes to a value of ∼ 90 % as σ ob decreases. For O 3 and N 2 O using the MLS precision values results in nearly the maximum WEP. For H 2 O, however, the MLS precision corresponds to a WEP value that is significantly smaller than the maximum. The fact that the WEP for all three tracers asymptotes to the same value suggests that the maximum WEP is limited by deficiencies in the DAS (e.g., TLM errors and imperfect modeling of the background errors), rather than limited by the geophysical variability of the tracers. To probe the theoretical limit of WEP, we also performed idealized 10-day simulations using unperturbed "perfect" observations and highly inflated background wind errors in order to artificially force larger wind increments. These experiments (not shown) resulted in maximum WEP of 98 % for O 3 and N 2 O and 95 % for H 2 O. This suggests that complete (100 % WEP) wind extraction is theoretically possible from these tracers, but the extraction in realistic experiments will be limited both by the DAS and by observation errors.
N 2 O provides the largest WEP values for a given relative observation error due to its large dynamic range, while the small dynamic range of H 2 O limits the wind extraction. As the error increases, the WEP values decrease monotonically. For all experiments shown, the effect of tracer assimilation is to improve the winds relative to the case of "no data" (indicated by black dotted line on Fig. 10). Even with imposed Gaussian random observation errors of 100 % of the global mean value, wind information is extracted in this 4D-Var system. Although the observation errors are large, the error is correctly specified in the assimilation system, so the data are correctly weighted. Therefore, any of these tracers has the potential to be useful for tracer wind extraction, even if they are poorly measured, as long as the data coverage is sufficient and the observation and background errors are modeled appropriately in the assimilation system. We caution, however, that these results are based on highly optimistic conditions. In a real world NWP system, tracer assimilation will be evaluated in terms of the incremental improvement to a wind analysis that is already fairly good. Model error, observation resolution, and observation biases might reduce or prevent any positive impact on winds.

Discussion
The use of trace gas data to benefit stratospheric winds in 4D-Var NWP systems is an attractive prospect, given the paucity of direct wind observations in the stratosphere. The success, however, depends on a number of factors such as the quality of the background tracer field and the observation precision, temporal sampling, and spatial coverage. Han and McNally (2010) reported that solar backscatter ultraviolet (SBUV) ozone assimilation could degrade the operational ECMWF 4D-Var wind analyses. In order to prevent erroneous wind increments, they stated that "the observation operator that links wind adjustments to changes in ozone concentration has been artificially cut." The erroneous increments were caused by the system attempting to correct for biases between ozone observations and the model background. Biases are one of the major obstacles that need to be overcome in order to extract wind from stratospheric tracers.
The main tracer that has been used in previous studies is O 3 , for which observations are readily available, and O 3 is also generally included as a prognostic variable in NWP forecast models. H 2 O is also included in NWP systems, so it provides an additional tracer for wind extraction. However, the results of this study show that the impact of H 2 O, when the observation errors are constrained by MLS precisions, on the dynamics in the middle stratosphere via the 4D-Var mechanisms is far less than for O 3 . N 2 O provides a third interesting possibility, which for these experiments rivals the potential benefits from O 3 . One advantage of N 2 O over O 3 is that the photochemistry of N 2 O is less important, given a long photochemical lifetime in the stratosphere. In this study, the tracers were assumed to be passive. Inclusion of O 3 and N 2 O errors in photochemical parameterizations could tip the scales towards N 2 O as being the more useful tracer for wind extraction at certain levels. However, it is likely that the impact of N 2 O would decrease with altitude in the upper stratosphere due to smaller mixing ratio and shorter photochemical lifetime.
This study modeled the wind-tracer coupling using a 4D-Var system, which has certain inherent limitations. In 4D-Var the coupling of tracer and dynamical fields (wind, temperature, and geopotential) occurs implicitly via the adjoint of the TLM, which propagates sensitivities of the cost function with respect to the tracer observations backwards in time to obtain sensitivities of the cost function to the model state at the beginning of the analysis window. The optimization problem involves perturbing the initial model state such that the resulting forecast minimizes the "distance" of the analysis with the background and with the observations, subject to error (and possibly other) constraints. In general, tracer-wind coupling is not specified explicitly in P b 0 . Realistic tracerwind coupling is more difficult to impose with a static P b 0 than the wind-temperature coupling. The latter takes advantage of dynamical balance approximations (Daley, 2001). No such balances are available between tracer and winds, although correlations between certain tracers and potential vorticity may prove useful (Li et al., 1998;Elbern et al., 2010).
One way to obtain realistic tracer-wind correlations in P b 0 is to use ensemble statistics to propagate information from one observed variable to the other ones. This approach avoids using the TLM and adjoint completely, relying on an ensemble of forecasts based on the full nonlinear forecast model. Milewski and Bourqui (2011) used an ensemble Kalman filter (EnKF) with a chemistry-climate model in a system that assimilated realistic stratospheric temperature and ozone observations. It was shown that assimilation of ozone profiles was able to significantly improve the wind analysis due to the cross-covariances. This suggests that a well-specified P b 0 may be used to constrain the wind field when assimilating tracer observations. Further work by Massart et al. (2012) examined using ensemble-based ozone P b 0 specification within a 4D-Var assimilation suite in order to improve ozone analyses. This "hybrid" approach is gaining momentum in the NWP community (e.g., Kuhl et al., 2013) but is relatively unexplored in the context of chemical data assimilation. While Massart et al. (2012) showed that ozone analyses benefitted from improved flow-dependent ozone P b 0 , their system was univariate (based on a chemistry and transport model driven by offline winds), and so the ozone assimilation did not provide feedback on the model dynamics. A full hybrid ensemble approach to the complete chemical and dynamical system may produce further improvements to tracer-wind extraction. Although the EnKF and hybrid-EnKF approaches for tracer-wind extraction are very promising, there are still many challenges that must be overcome, not least of which is the computational cost. Therefore, many operational NWP centers (including the US Navy's operational Fleet Numerical Meteorology and Oceanography Center) are still currently using 4D-Var assimilation techniques.

Conclusions
A 4D-Var data assimilation system based on a global shallow water model coupled to a tracer advection equation was used to examine the tracer-wind coupling problem in a highly idealized setting. While a number of additional tracers (methane, chlorofluorocarbons, etc.) could have been tested, we focused on O 3 , N 2 O, and H 2 O, which are readily available from instruments such as MLS. Using tracer data based on a nature run, we assimilated hourly observations on the full Gaussian grid. In experiments with very small observation error (BEST cases), the three tracers are able to constrain the model winds, height, and tracer, with wind extraction potential of 90 % after 10 days, wind errors of ∼ 0.3 m s −1 , and height errors of ∼ 13 m, given extremely small observation errors of 0.0419 ppmv for O 3 , 0.00274 ppmv for N 2 O, and 0.0157 ppmv for H 2 O. For the BEST cases, each tracer was able to provide similar dynamical information. However, when the tracer observations have errors that match the precision of MLS, O 3 provides the most dynamical information, closely followed by N 2 O. Due to weaker gradients relative to the MLS precision, H 2 O provided the least impact on the winds. As expected, wind extraction is more difficult with larger observation error, although the overall impact on the winds is positive even with random errors up to 100 % of the global mean, as long as the observation error standard deviation is specified consistently with the imposed errors. These results show that by taking advantage of the tracer-wind coupling mechanism, 4D-Var assimilation of global tracer observations is able to very accurately constrain the entire state vector of the shallow water model (winds, height, and tracer) under idealized conditions. In a realistic NWP system with a full suite of observations, errors in the forecast model and tracer observation biases and sampling limitations will likely reduce the impact of tracers on the dynamics.