Abstract

A methodology is presented for objectively optimizing nonorographic gravity wave source parameters to minimize forecast error for target regions and forecast lead times. In this study, we employ a high-altitude version of the Navy Global Environmental Model (NAVGEM-HA) to ascertain the forcing needed to minimize hindcast errors in the equatorial lower stratospheric zonal-mean zonal winds in order to improve forecasts of the quasi-biennial oscillation (QBO) over seasonal time scales. Because subgrid-scale wave effects play a large role in driving the QBO, this method leverages the nonorographic gravity wave drag (GWD) parameterization scheme to provide the necessary forcing. To better constrain the GWD source parameters, we utilize ensembles of NAVGEM-HA hindcasts over the 2014–16 period with perturbed source parameters and develop a cost function to minimize errors in the equatorial lower stratosphere compared to analysis. Thus, we may determine the set of GWD source parameters that yields a forecast state that most closely agrees with observed QBO winds over each optimization time interval. Results show that the source momentum flux and phase speed spectrum width are the most important parameters. The seasonal evolution of optimal parameter value, specifically a robust semiannual periodicity in the source strength, is also revealed. Changes in optimal source parameters with increasing forecast lead time are seen, as the GWD parameterization takes on a more active role as QBO driver at longer forecast lengths. Implementation of a semiannually varying source function at the equator provides RMS error improvement in QBO winds over the default constant value.

1. Introduction

It has long been understood that internal gravity waves play an integral role in the large-scale circulation of the middle atmosphere. They are emitted from a variety of tropospheric sources, both orographic and nonorographic. Some propagate into the stratosphere and mesosphere, where they break and deposit some or all of their momentum locally. Since this process occurs quasi continuously around the globe, it sustains planetary-scale gravity wave drag (GWD) forces that drive the climate and meteorology of the middle atmosphere across a broad range of spatial and temporal scales (Fritts and Alexander 2003).

While both orographic and nonorographic gravity waves are critical drivers of the extratropical stratosphere and mesosphere (e.g., Holton 1983; McLandress 1998; Hitchman et al. 1989; Alexander and Rosenlof 2003; Scinocca et al. 2008), the dominant quasi-biennial and semiannual circulations of the tropical stratosphere and mesosphere are predominantly impacted by nonorographic GWD (e.g., Dunkerton 1982, 1997; Baldwin et al. 2001; Peña-Ortiz et al. 2010). Since finite computational resources force current global numerical weather prediction (NWP) and climate models to operate at space–time resolutions that do not resolve short-wavelength components of the atmospheric wave spectrum, accurate parameterizations of subgrid-scale nonorographic GWD remain an indispensable component of these models.

Most nonorographic gravity wave drag parameterization schemes begin by specifying the source-level wave features based on observational estimates of momentum flux and spectral properties in the real atmosphere, which are used to build a parameterized source-level wave spectrum whose propagation upward and subsequent dissipation are modeled by a chosen dynamical framework for parameterizing these processes. Typically, the global atmospheric wave spectrum, in particular the short wavelength regime represented by the GWD parameterization, has been characterized using decades of radiosonde and satellite observations (Allen and Vincent 1995; Tsuda et al. 2000; Alexander et al. 2010; Zhang et al. 2012). In contrast to curating climatological nonorographic gravity wave spectra through offline tuning, some studies have proposed and implemented schemes where source characteristics are derived directly from parameterized subgrid-scale convective activity within the model, which is a more internally consistent approach based on our physical understanding of gravity wave sources (e.g., Beres et al. 2004; Richter et al. 2010; Schirber et al. 2014). However, varying implementations of nonorographic GWD parameterization in models and differences in modeled large-scale flow patterns mean that the parameterized wave spectrum that provides the highest NWP or climate skill may differ among models (e.g., Bushell et al. 2015; Serva et al. 2018).

A persistent practical difficulty with all parameterizations of nonorographic GWD is the need for ad hoc adjustments (tuning) once implemented in a model, a problem further complicated as models enter a new “gray zone” era where both gravity waves and convection are partially resolved and partially parameterized (e.g., Vosper et al. 2016; Chen et al. 2018). As horizontal and vertical resolutions are routinely changed in any given model from application to application, the partitioning between parameterized and resolved GWD also changes, necessitating careful retuning of the parameterized GWD component. Since the amount of resolved gravity wave momentum flux (and hence the corresponding amount of GWD requiring parameterization) also proves sensitive to many other model characteristics besides grid resolution (e.g., the levels and types of numerical diffusion; see Skamarock 2004; Shutts and Vosper 2011), objective procedures for retuning parameterized GWD in response to changes in model properties have proven elusive. To the limited extent that offline GWD tuning methods for the middle atmosphere have been discussed, they focus mainly on ad hoc “trial and error” approaches in which GWD parameters are exhaustively modified to most closely reproduce some subset of observational states (see, e.g., Kim et al. 2003; Scinocca et al. 2008; Eckermann et al. 2009; Long et al. 2014).

A more promising approach to emerge recently uses objective data assimilation methods to identify subgrid-scale drag deficits in model forecasts, which are then minimized using a cost function approach that depends on key parameters of the model’s GWD parameterization (e.g., Pulido et al. 2012; Tandeo et al. 2015). We explore this type of objective approach here for the specific problem of tuning parameterized nonorographic GWD driving of the tropical stratosphere within a model.

Nonorographic GWD parameterization is an important component of an accurate simulation or forecast of the tropical quasi-biennial oscillation (QBO), which is known to be primarily driven by vertically propagating waves (Lindzen and Holton 1968; Holton and Lindzen 1972; Baldwin et al. 2001) with a significant contribution from small-scale nonorographic gravity waves (Alexander and Holton 1997; Dunkerton 1997). The QBO is also a potent source of predictability in the middle atmosphere owing to relatively long tropical stratospheric memory (Scaife et al. 2014). Many studies have characterized the impact of the QBO on forecast skill in regions across the globe (Holton and Tan 1980; Hitchman and Huesmann 2009; Garfinkel and Hartmann 2011; Garfinkel et al. 2018). Since accurate simulation of the QBO is important for extended-range prediction, tuning a model’s nonorographic GWD parameterization specifically for the QBO is an important step toward this goal. The present work is timely given new coordinated initiatives to improve QBO modeling and prediction, and in particular the role of parameterized nonorographic GWD (Butchart et al. 2018).

We propose in this study a methodology for determining the optimal source parameters of a nonorographic GWD parameterization, in the sense that forecast error is minimized. In particular, we are interested in improving seasonal forecasts of the stratospheric QBO winds, but the methodology may be adapted to regions of any size and location globally. Details of the working forecast model, its nonorographic GWD parameterization, and the formulation of the optimization technique are given in section 2. In section 3 we describe our experimental design, including model improvements and specification of the methodology for accurate QBO simulation. Results of the optimization are shown in section 4, and improvements in forecast skill are discussed in section 5. A summary of the study’s findings follows in section 6.

2. Model and methodology

a. NAVGEM-HA model description

In this study, analysis and forecast winds were provided by a high-altitude (HA) configuration of the Navy Global Environmental Model (NAVGEM-HA; see, e.g., McCormack et al. 2017). As the Navy’s operational global NWP system, NAVGEM couples a global model issuing forecasts with a data assimilation system (DAS) issuing meteorological analyses that serve both as atmospheric initial conditions and verification fields for the forecasts (Hogan et al. 2014).

For this work, we operated NAVGEM-HA with triangular spectral truncation at total wavenumber 119 (T119, corresponding to gridbox dimensions of ~1° on the quadratic Gaussian grid), with 74 (for the analysis) or 108 (for the forecasts) hybrid σp model levels. As shown in Fig. 1, the L108 configuration provides improved vertical resolution within the 100–1-hPa region of interest in this work, relative to other layer profiles used in recent NAVGEM studies (e.g., Kuhl et al. 2013; Barton and McCormack 2017; Eckermann et al. 2018). Most notably for this study, L108 provides the requisite stratospheric vertical resolution to sustain a prognostic QBO in the model (e.g., Giorgetta et al. 2006; McCormack et al. 2015; Anstey et al. 2016).

Fig. 1.

NAVGEM pressure–height layer thicknesses ΔZk vs pressure altitude Zk for L108 (red) used here, operational L60 (black) (Kuhl et al. 2013), L74 HA (green) (Eckermann et al. 2018), and L80 (blue) (Barton and McCormack 2017). The gray band shows the altitude region focused on in the present study.

Fig. 1.

NAVGEM pressure–height layer thicknesses ΔZk vs pressure altitude Zk for L108 (red) used here, operational L60 (black) (Kuhl et al. 2013), L74 HA (green) (Eckermann et al. 2018), and L80 (blue) (Barton and McCormack 2017). The gray band shows the altitude region focused on in the present study.

The T119L74 NAVGEM-HA configuration was run in an historical reanalysis mode to generate a contiguous set of meteorological reanalyses spanning a 3-yr period from 2014 to 2016. These NAVGEM-HA analyses were produced using a four-dimensional variational (4DVAR) data assimilation (DA) algorithm, with background errors specified as linear combinations of static covariances and flow covariances from an 80-member forecast ensemble at the inner loop (T47) resolution (so-called hybrid 4DVAR; see Kuhl et al. 2013). In the 100–1-hPa range of interest, the system assimilated data from a dense network of heterogeneous satellite and suborbital observations (see, e.g., Figs. 6–8 of Eckermann et al. 2018). Following McCormack et al. (2017), we combined the resulting 6-h analyses with the 3-h outer loop forecast from each update cycle to produce global fields every 3 hours, to reduce tidal aliasing at upper levels.

These reanalysis fields are used both as initial conditions and verification for extended-range T119L108 forecasts used here to optimize parameterized nonorographic GWD. The NAVGEM forecast model used here is global, semi-implicit and semi-Lagrangian, and incorporates a suite of parameterizations tuned for tropospheric and stratospheric NWP (Hogan et al. 2014). The HA configuration is supplemented with dynamical modifications and additional physical parameterizations for mesospheric NWP, including parameterized nonorographic GWD and exothermic chemical heating (e.g., McCormack et al. 2017).

b. NAVGEM-HA GWD parameterization

Because the specifics of the nonorographic GWD parameterization scheme are relevant to our optimization approach described in later sections, we briefly review its formulation in NAVGEM-HA. The basic approach, illustrated schematically in Fig. 2, follows closely that of the version 3.0 Whole Atmosphere Community Climate Model (WACCM; Garcia et al. 2007). At each horizontal grid point, a spectrum of waves is simulated, with each wave assigned a unique ground-based horizontal (zonal) phase speed cj in the range usrc ± cr, where usrc is the zonal wind at the launch level (model pressure psrc) and cr is called the phase speed range. By design, the scheme does not parameterize waves with source-level intrinsic phase speeds |cjusrc| greater than |cr|. The source-level momentum flux of wave j is a Gaussian function of intrinsic phase speed:

 
τj=τ×exp[(cjusrc)2c02],
(1)
 
τ=τsrcF(φ,t),

where τ is the prescribed peak source flux, specified using a base value τsrc modulated by a latitudinally and seasonally varying component F(φ, t), and c0 is the phase speed width. The Gaussian is truncated at cj = usrc ± cr (see Fig. 2, bottom panel). By centering the Gaussian at usrc, wave momentum is distributed symmetrically at the source level among waves with positive (eastward) and negative (westward) intrinsic phase speeds.

Fig. 2.

Schematic illustration of the nonorographic GWD scheme in NAVGEM-HA. A zonal wind profile for an equatorial grid point is shown, and two waves with phase speeds cj = ±30 m s−1 are launched (green curves) at psrc ~ 500 hPa. The bottom plot shows the Gaussian variation of source-level momentum flux with phase speed given source-level zonal wind usrc ~ 4 m s−1, τsrc ~ 4.5 mPa, c0 ~ 30 m s−1, and cr ~ 50 m s−1. The left and right panels demonstrate how each wave is attenuated at levels where the wave momentum flux (green) exceeds the local saturation limit calculated from Eq. (2) (blue). The left panel uses signed flux values for illustration only. Corresponding zonal wind acceleration [Eq. (3)] due to wave breaking is plotted in red.

Fig. 2.

Schematic illustration of the nonorographic GWD scheme in NAVGEM-HA. A zonal wind profile for an equatorial grid point is shown, and two waves with phase speeds cj = ±30 m s−1 are launched (green curves) at psrc ~ 500 hPa. The bottom plot shows the Gaussian variation of source-level momentum flux with phase speed given source-level zonal wind usrc ~ 4 m s−1, τsrc ~ 4.5 mPa, c0 ~ 30 m s−1, and cr ~ 50 m s−1. The left and right panels demonstrate how each wave is attenuated at levels where the wave momentum flux (green) exceeds the local saturation limit calculated from Eq. (2) (blue). The left panel uses signed flux values for illustration only. Corresponding zonal wind acceleration [Eq. (3)] due to wave breaking is plotted in red.

The upward vertical propagation of waves from this source and flux deposition at upper levels is parameterized, following Lindzen (1981), using hydrostatic wave equations that ignore any lateral group components. At each model level in the column above the launch height, the momentum flux saturation threshold for wave j is calculated as

 
τsat,j=|kρ(ucj)3Frc2N2|,
(2)

where k is zonal wavenumber, ρ is density, u is the background zonal wind at that level, Frc is the critical Froude number for wave breaking (a constant of order unity), and N is the buoyancy frequency at that level. If the momentum flux carried by a wave exceeds this local saturation limit, the excess is deposited into the background flow. The resulting momentum flux divergence provides a zonal mean-flow acceleration:

 
aj=gϵτjp,
(3)

where g is gravitational acceleration and ϵ is a parameterized efficiency (or intermittency) of wave breaking. The acceleration aj is given the same sign as ucj. Note that a parameterized wave is totally absorbed at its critical level where cj = u because τsat,j = 0. In the left and right panels of Fig. 2, we show how the saturation limit may vary with height and the wave’s phase speed, especially in relation to the background zonal wind profile, and how the momentum carried vertically by the parameterized waves is deposited.

The total GWD tendency is obtained by summing the aj from all parameterized waves j in Eq. (3). Since a large number of waves is required to discretize the flux Gaussian in Eq. (1) accurately, this GWD algorithm can be computationally expensive when implemented in models. Eckermann (2011) modified this scheme into a stochastic analog to improve both the computational efficiency and ensemble spread of GWD within models, which facilitated its transition into NAVGEM for use in operational NWP. This stochastic version, used in NAVGEM-HA and hence in this work, treats only a single wave at each grid point, with the phase speed randomly assigned in the range usrc ± cr. The efficiency ϵ must be upscaled to account for the difference in total source momentum flux between the multiwave and single-wave schemes, but otherwise no other changes are necessary since, as shown by Eckermann (2011), the same time-mean flux and mean-flow acceleration are produced, leading to similar time-mean circulation responses in models. A range of subsequent studies have pursued similar stochastic approaches to parameterizing nonorographic GWD and have reported similarly acceptable performance in modeling both tropical and extratropical middle atmospheric circulations (e.g., Lott and Guez 2013; McCormack et al. 2015; de la Cámara et al. 2016; Garcia et al. 2017; Serva et al. 2018).

While the basic approach of the GWD parameterizations used in NAVGEM-HA and WACCM 3.0 (Garcia et al. 2007) are similar, the implementation of the scheme in NAVGEM-HA differs from the WACCM 3.0 implementation in a number of ways. Effects of radiative damping and eddy diffusive mixing are turned off, and all parameterized momentum that reaches the upper sponge layer of the model is deposited there so that circulation responses due to downward control may be realized (Haynes et al. 1991; Garcia and Boville 1994). The GWD-induced heating rates are calculated differently [see Eq. (7) of Eckermann 2011]. Although the temperature tendency is applied by default in both WACCM 3.0 and NAVGEM-HA, for this study we did not apply this tendency, since this term is generally important only in the mesosphere and higher.

The standard NAVGEM-HA GWD parameterization is tuned so that momentum is deposited mostly in the extratropical mesosphere. Motivation for the present study arose from the realization that our T119L108 model does not resolve the full spectrum of convectively generated tropical waves needed to drive the QBO (e.g., Krismer et al. 2015), that nonorographic GWD parameterized using the default configuration does not generate and sustain a realistic QBO in free-running model simulations, and that an objective retuning of the source parameters (e.g., McCormack et al. 2015) would be required inter alia to accurately forecast the QBO in equatorial stratospheric winds on seasonal time scales. To this end, we have selected five parameters for optimization in this study: “source strength” τsrc (default 4.5 mPa), launch height psrc (default 500 hPa), phase speed Gaussian width c0 (default 30 m s−1), phase speed range cr (default 80 m s−1), and zonal wavenumber k (default wavelength 2π/k = 400/a for radius of Earth a). These parameters are assigned globally, while the source flux varies around the nominal τsrc = 4.5 mPa value with latitude and season according to a prescribed F(φ, t) in Eq. (1) that varies with season in the extratropics but is seasonally invariant in the tropics (see Fig. 1b of Eckermann 2011).

The specific influence of each source parameter on the parameterized acceleration aj varies and so can be used to tune circulation responses in the model (see, e.g., section 3 of Eckermann et al. 2009). Changes in the source strength τsrc produce a uniform relative change in source flux across the parameterized spectrum in Eq. (1), meaning each wave j will see an equivalent scaling of its source-level momentum flux leading to more GWD throughout model upper layers. However, changes to the phase speed width do not affect waves across the entire spectrum equally; in fact, only medium-speed (|cjusrc| ~ 30 m s−1) waves will see a significant absolute change in source momentum flux, while waves with phase speeds nearer usrc or in the tails of the distribution receive little to no absolute change in source momentum flux. In general, raising the source flux by any means in isolation will produce larger parameterized accelerations and will tend to move the onset of wave–mean flow interaction to lower levels. The phase speed range eliminates any waves faster than a desired intrinsic phase speed limit and therefore has little effect when this value is larger than typical wind speeds in the model layers of interest, because the waves at the ends of the spectrum are too fast to interact with the background flow locally. Because wavenumber k acts more or less as a simple scaling coefficient for the saturation flux profile via Eq. (2), a higher value will tend to move breaking altitudes up, meaning, for example, that more flux is deposited nearer to a critical level, leading to more sharply peaked aj profiles. Last, the launch height controls the depth of the filtering layer below the tropopause, which can lead to significantly different GWD profiles and middle-atmosphere responses (e.g., Manzini and McFarlane 1998).

The goal of this study is to optimize the GWD source parameters to produce realistic forecasts of the QBO. In the standard NAVGEM-HA implementation, source parameters are assigned one value for the entire global domain, so changing their values may have a more drastic effect on the overall circulation of the forecast model than desired, which could indirectly affect the QBO and negatively impact the fidelity of the optimization technique. We have therefore modified the nonorographic GWD parameterization scheme in NAVGEM-HA so that the five parameters selected for optimization (τsrc, psrc, c0, cr, and k) are varied independently within latitude bands. The parameters are now specified at a set of “anchor latitudes,” which we have chosen to be at 90°, 60°, 40°, and 20°S; 0°; and 20°, 40°, 60°, and 90°N, and vary linearly between adjacent latitudes. In the equatorial band 5°S–5°N, the values are kept equal to the specified equatorial anchor point value rather than being linearly interpolated. Unless noted otherwise, for ensembles of QBO hindcasts analyzed in this study the parameters are given their default value at each anchor latitude, and we apply perturbations in the equatorial band only.

c. GWD parameter optimization

Here we develop a general methodology by which we can determine the “best” (optimal) values of a subset of GWD source parameters that minimize the model’s hindcast error, relative to analysis, in equatorial stratospheric zonal-mean zonal wind u¯. The simplest tuning method, often called ad hoc tuning, is to run a large number (i.e., an ensemble) of hindcasts with varying source parameters and select the GWD tuning of the member with the smallest error. However, with a large set of tuning parameters, the probability of sampling the optimal parameter space is small, in addition to being inefficient.

The problem can be simplified if the hindcast error at all altitudes of interest can be linearly related to changes in a GWD source parameter over some range of perturbations. Assuming a change in the u¯ profile over the hindcast period Δt of Δu¯, we can then express the model response as a simple linear matrix M governing the perturbation equation:

 
Mx=Δu¯=Δu¯Δu¯def,
(4)

where x′ is the vector of perturbations of the GWD source parameters from their default values and Δu¯def is the change in u¯ over the hindcast period Δt for default values of the GWD parameters. Inclusion of this term satisfies the trivial case x′ = 0. Because the hindcasts are launched from the same set of initial conditions at time t0, Δu¯ simplifies to (u¯u¯def)|t=t0+Δt=u¯, which is just the difference in hindcast wind profiles at the final time t. Thus, M provides the response of the hindcast zonal-mean zonal wind profile at lead time Δt to changes in the GWD source parameters. Each column in M quantifies the linearized Δu¯ response to perturbations in one of the GWD source parameters, while each row in M quantifies the response within a different vertical model layer.

To minimize the total error between the hindcast and our verifying analysis (hereafter the AF error) we first replace u¯|t=t0+Δt above with the state we require, u¯A, the QBO wind profile from our verifying analysis at the final time t. Since the linear matrix equation in Eq. (4) is overdetermined, with more analyzed wind values u¯A on model layers (rows) to reproduce than source parameters x′ available to vary (columns), an exact x′ solution does not exist. Thus, we estimate a “best” x′ from Eq. (4) using least squares minimization of the norm Mx(u¯Au¯def)|t=t0+Δt.

One implicit assumption of such an approach is that all of the AF error is due to forecast error from the nonorographic GWD parameterization scheme. While this is never the case in practice (see section 3a), with respect to the QBO it is likely that the majority of the wind error comes from insufficient wave driving by the parameterized nonorographic gravity waves (Scaife et al. 2000; Giorgetta et al. 2002; Geller et al. 2016). As a consequence, this methodology will minimize the model error from any source so long as that error projects strongly onto the model response due to changes in nonorographic GWD source parameters.

Since Δu¯ is the result of integration of the full nonlinear forecast model, M must contain information about the change in zonal wind due to the combined effects of every model component, not just the GWD parameterization. While this makes it infeasible to obtain M precisely, we can estimate it as follows. First, we note that M as defined in Eq. (4) is conceptually the same tangent-linear form of the forecast model utilized in the 4DVAR DA algorithm [e.g., Eq. (6) of Eckermann et al. 2018], but here the tangent-linear model (TLM) is one linearized around a background forecast trajectory defined by the default source-level parameters of the GWD parameterization. Thus, we can populate all the elements of the TLM matrix M by first performing an ensemble of hindcasts in which the GWD source parameters are perturbed symmetrically about the default state. Matrix M is the matrix that best satisfies (in a least squares sense) the following system of equations:

 
M[x1xE]X=[u¯1u¯E]U,
(5)

where each column in X is the vector of perturbations of the GWD source parameters in the selected subset for ensemble member e (xe), and the corresponding column in U is the deviation of the hindcast profile of ensemble member e (u¯e) from the default hindcast (u¯def) at time t = t0 + Δt. We solve Eq. (5) using multiple linear regression to fit the linear model M in the multidimensional parameter space.

To gauge the validity of the linearity assumption inherent to M for a given GWD source parameter subset, we correlate the perturbation of the actual ensemble member hindcast profile, u¯e, with the perturbation profile predicted using the linear model M and the ensemble member parameters xe [i.e., u¯pred=Mxe as in Eq. (4)]. The squared correlation coefficient (R2) between u¯e and u¯pred will be used to assess the accuracy of the linearization approximation, with high R2 values denoting strong linearity. Low R2 values may indicate that interaction terms, which we assumed did not exist when we wrote Eq. (4), contribute significantly to u¯. Poor prediction skill in the region of interest may also promote low R2 by introducing significant hindcast error that is unrelated to the GWD source parameters. Examples and validation of the approach will be provided in the next section.

3. Experimental setup

a. Model adjustments for QBO simulation

The change in zonal mean wind over the forecast period Δu¯ can be directly influenced by other aspects of the model apart from parameterized GWD, such as the dynamical core (Yao and Jablonowski 2015), horizontal and vertical resolution (Hamilton et al. 1999; Giorgetta et al. 2006; Anstey et al. 2016), convection parameterizations (Horinouchi et al. 2003; Yang et al. 2011), and vertical diffusion (McCormack et al. 2015). Here, we describe how these components of the forecast model are configured to support a realistic QBO simulation in NAVGEM-HA hindcasts that produce the background forecast trajectory Δu¯def around which our governing GWD matrix equation in Eq. (4) is linearized.

An ensemble of T119L108 hindcasts is generated by vertically interpolating our T119L74 analysis fields linearly with respect to ln(p) to L108 (see Fig. 1) for use as initial conditions. There is no extrapolation required in the vertical regridding because L74 has a higher model top than L108, and we expect that noise introduced by initializing the model with the interpolated fields will be small. Over the hindcast run, sea surface temperature and sea ice are updated every 12 forecast hours using archived operational analyses provided by Navy Fleet Numerical Meteorology and Oceanography Center. The tendency due to vertical diffusion is reduced by two orders of magnitude in the stratosphere alone to preserve vertically narrow (<1 km) QBO shear zones [see McCormack et al. (2015) for details of the NAVGEM-HA vertical diffusion scheme and its impact on simulating the QBO]. Last, we have replaced the default simplified Arakawa–Schubert cumulus convection scheme (SAS; see Hogan et al. 2014) with a modified Kain–Fritsch scheme (Ridout et al. 2005), which produces a more realistic resolved wave spectrum in the tropical upper troposphere and lower stratosphere over forecast periods longer than 5 days (Janiga et al. 2018).

b. Single-parameter experiments

Figures 37 detail results from the ensemble of equatorial u¯ hindcasts, where only the indicated GWD parameter has been varied among ensemble members while the other four parameters are kept at their default value. This is consistent with single columns of the tangent-linear matrix M that in turn facilitate direct comparisons with the linear matrix predictor [Eq. (4)]. Here we may posit linearity of the single-parameter responses before exploring the linearity of M in the multidimensional parameter space. Hindcast ensembles are initialized at 0000 UTC 1 March 2014 and run for 30 days, and so are compared in these figures to the verifying analysis valid at 0000 UTC 31 March 2014 (black curves). The left panels in each figure show the zero-wind line as a function of pressure and time for the analysis (black) and ensemble members (color-coded according to parameter value), while the middle panels show the vertical profiles of u¯ at the 30-day lead time. The scatterplots in the right panels correlate the perturbation of the actual response from hindcast ensemble member e (u¯e) with the predicted response u¯pred=Mxe as in Eq. (4) over all analysis levels between 70 and 1 hPa (16 levels).

Fig. 3.

Ensemble spread of 30-day hindcasts with varying phase speed width c0 within the range ±20 m s−1 from the default value of 30 m s−1 (see color bar). All hindcasts were initialized at 0000 UTC 1 Mar 2014 and run out for 30 days. (a) Zero-wind line as a function of time and pressure for the NAVGEM-HA analysis u¯A (black) and ensemble hindcast members (colored according to c0 value). (b) Final vertical profiles of u¯ and u¯A valid at 0000 UTC 31 Mar 2014 (30-day lead time). (c) Scatterplot correlating u¯ predicted from Eq. (4) using M and xe and the actual hindcast perturbation profile u¯e at the final time over all analysis levels between 70 and 1 hPa (16 levels).

Fig. 3.

Ensemble spread of 30-day hindcasts with varying phase speed width c0 within the range ±20 m s−1 from the default value of 30 m s−1 (see color bar). All hindcasts were initialized at 0000 UTC 1 Mar 2014 and run out for 30 days. (a) Zero-wind line as a function of time and pressure for the NAVGEM-HA analysis u¯A (black) and ensemble hindcast members (colored according to c0 value). (b) Final vertical profiles of u¯ and u¯A valid at 0000 UTC 31 Mar 2014 (30-day lead time). (c) Scatterplot correlating u¯ predicted from Eq. (4) using M and xe and the actual hindcast perturbation profile u¯e at the final time over all analysis levels between 70 and 1 hPa (16 levels).

Fig. 4.

As in Fig. 3, but for source strength τsrc varied within the range ±2 mPa around the default value of 4.5 mPa.

Fig. 4.

As in Fig. 3, but for source strength τsrc varied within the range ±2 mPa around the default value of 4.5 mPa.

Fig. 5.

As in Fig. 3, but for phase speed range cr varied within the range 10–90 m s−1.

Fig. 5.

As in Fig. 3, but for phase speed range cr varied within the range 10–90 m s−1.

Fig. 6.

As in Fig. 3, but for launch height psrc varied within the range 100–500 hPa.

Fig. 6.

As in Fig. 3, but for launch height psrc varied within the range 100–500 hPa.

Fig. 7.

As in Fig. 3, but for wavenumber k varied within the range 100/a–900/a.

Fig. 7.

As in Fig. 3, but for wavenumber k varied within the range 100/a–900/a.

Phase speed width (Fig. 3) and source strength (Fig. 4) both exemplify a good predictor for the purpose of using our linear matrix model [Eq. (4)] to relate changes in a single predictor x′ to a corresponding u¯. Here we define a good predictor as a source parameter with a large ensemble spread that exhibits a zonal wind response that varies linearly with the parameter value. A parameter with these characteristics ensures that our assumption of linearity required for Eq. (4) to be accurate holds, and thus that the optimization algorithm will return a realistic and reasonable optimal value for this parameter. Figures 3 and 4 display a smooth, monotonic change in wind with respect to parameter value. Note that below 70 hPa, it is evident from the lack of monotonicity that the linearity assumption is invalid, thus we will take 70 hPa as the domain’s lower bound for optimization fits. For this particular case, our least squares solution to the matrix Eq. (4) yields optimal values of c0 and τsrc of about 24 m s−1 and 3.0 mPa, respectively. The R2 values are calculated to be 0.941 and 0.987, which again verifies that a single-predictor model using c0 or τsrc is able to accurately reproduce hindcast profiles of u¯ using the perturbed parameter value and our tangent-linear matrix M. Note also that changes in c0 and τsrc affect the descent rate of the QBO more strongly than its amplitude in 30-day hindcasts, with smaller parameter values corresponding to slower descent.

The total source momentum flux, which is the integral of the Gaussian in Eq. (1), is proportional to the product τsrcc0. Equivalent changes to τsrc and c0 therefore change the available momentum flux at the launch height equally. However, some of the difference in total source momentum flux created by changing τsrc, which uniformly increases or decreases the source flux of all waves by the same multiplicative constant, is not available to use in the equatorial stratosphere. Power added or removed from the fastest or slowest waves will have little to no effect on the winds in our domain of interest because those waves are dissipated above the 1-hPa level or below the 100-hPa level. On the other hand, recall that the effect of changes to c0 scales relative to the parameterized wave’s intrinsic phase speed. Thus, we expect that changes in c0 (Fig. 3) will produce a larger ensemble spread than τsrc (Fig. 4) because of a larger change in available momentum flux applied to those parameterized waves most relevant for the QBO.

Phase speed range cr behaves uniquely among the five source parameters (Fig. 5), in that it appears to exhibit linearity only within a small window between cr = 20 and 40 m s−1, while outside of this range the ensemble member curves reproduce roughly the same wind profiles found using the bounding cr values. Consider the case of cr < 20 m s−1, meaning only waves with intrinsic phase speeds less than 20 m s−1 are parameterized. In this case, it is likely that most of the parameterized waves are dissipated in the upper troposphere and lower stratosphere before reaching QBO altitudes, so that changes to cr here have no effect on the hindcast QBO winds. Likewise, introducing waves into the parameterization with phase speeds greater than about 40 m s−1 will have no effect on the hindcast QBO winds because these waves are generally too fast to interact with typical background winds in the equatorial lower stratosphere. Taking these results into account, we believe that cr is not a good choice for prediction of QBO winds by itself, although we note that it is indeed an important parameter elsewhere in the atmosphere and may be relevant for GWD tuning to simulate the stratospheric and mesospheric semiannual oscillations (SAO; see Fig. 5b), where nonorographic GWD from faster waves is known to play an important role (e.g., Garcia et al. 1997).

Launch height psrc (Fig. 6) and horizontal wavenumber k (Fig. 7) are also not good predictors mainly because they exhibit a relatively small ensemble spread. It is clear from both figures that the effects of psrc and k on the hindcast u¯ profiles are too minimal for a realistic value of either parameter to provide a measurable reduction in AF error. Note that this conclusion refers only to errors in u¯ in the equatorial stratosphere for hindcasts of subseasonal length and does not preclude either parameter from having a significant impact in other regions of the atmosphere or for longer lead times. For example, decreasing k in the tropics may be necessary for successful multiyear simulation of the QBO (Xue et al. 2012; McCormack et al. 2015; Yu et al. 2017), especially since the importance of inertia–gravity waves with wavelengths greater than 100 km for the evolution of the QBO has been shown (Kawatani et al. 2010; Holt et al. 2016).

c. Multiparameter experiment

To further explore the optimal GWD source parameter settings in the multidimensional parameter space composed of τsrc, psrc, c0, cr, and k, we utilize a random “cloud” method in which we run a 60-member ensemble of hindcasts, and for each ensemble member apply a random perturbation to every GWD source parameter in the selected subset. This differs from the standard TLM approach in Eq. (4) in which the optimized source parameter vector is found by running single-parameter optimizations (as in Figs. 37) iteratively until a threshold cost function minimum is reached. We utilize this method here because it allows us to quickly explore possible nonlinear interactions among predictors. If the nonlinear interaction terms have a significant impact on u¯, we would find low predictive capability of the simple linear model M described in the previous section [Eq. (4)] because it does not account for these interaction terms by design. We determine the linear model M that best fits perturbations in the multidimensional parameter space to the hindcast equatorial zonal-mean zonal wind response across all model levels between 70 and 1 hPa by solving Eq. (5) using multiple linear regression. The source parameter vector x′ composed of all five selected parameters is then optimized following the procedure outlined in section 2c.

The results of this experiment are shown in Fig. 8. Each blue curve in Figs. 8a–e is the profile of response of the hindcast equatorial zonal mean zonal wind (i.e., Δu¯) to a normalized unit change in the indicated source parameter. In Fig. 8f, the response predicted by the model M to the optimized GWD parameter vector x* [i.e., Δu¯=u¯optimalu¯def=Mx* as in Eq. (4)] is shown in red, and the validating analysis Δu¯ is in black. We find good agreement between these two profiles, which would indicate that the fitted linear model M is providing an accurate fit to the hindcast results. However, there is high correlation among the responses of all source parameters (multicollinearity) throughout the depth of the domain, notwithstanding the psrc response at heights above about the 7-hPa level. Note that multicollinearity does not preclude a model from accurately predicting the total response (cf. Fig. 8f), but it does inhibit our ability to separate the responses due to individual predictors. As such, we are unable to simultaneously estimate the optimal value for each of the individual parameters with a high degree of confidence. In the scatterplot (Fig. 8g), we find that the correlation between the predicted perturbation profile u¯pred=Mxe and the actual hindcast ensemble member perturbation profile u¯e has an R2 value of 0.554, indicating only moderate confidence that the model could closely reproduce the hindcast wind profile u¯e of ensemble member e when estimated using only M and xe as in Eq. (4). Since no pair of predictors is uncorrelated and nonlinear terms are likely to affect u¯ significantly, a model using only a single parameter (i.e., Figs. 37) is the best choice for the optimization procedure when applied to equatorial stratospheric u¯.

Fig. 8.

Response of 30-day hindcast u¯ wind profile initialized at 0000 UTC 1 Mar 2014 to a (normalized) unit change in (a) source strength, (b) launch height, (c) phase speed width, (d) phase speed range, and (e) wavenumber. (f) Change in analysis u¯ over the 30-day hindcast period (black) and total response of 30-day hindcast wind profile with optimized gravity wave drag source parameters (red). (g) Scatterplot correlating u¯ predicted from Eq. (4) using M and xe and the actual hindcast perturbation profile u¯e at the final time over all analysis levels between 70 and 1 hPa (16 levels).

Fig. 8.

Response of 30-day hindcast u¯ wind profile initialized at 0000 UTC 1 Mar 2014 to a (normalized) unit change in (a) source strength, (b) launch height, (c) phase speed width, (d) phase speed range, and (e) wavenumber. (f) Change in analysis u¯ over the 30-day hindcast period (black) and total response of 30-day hindcast wind profile with optimized gravity wave drag source parameters (red). (g) Scatterplot correlating u¯ predicted from Eq. (4) using M and xe and the actual hindcast perturbation profile u¯e at the final time over all analysis levels between 70 and 1 hPa (16 levels).

d. Methodology validation

Large R2 values calculated using our linear matrix model [Eq. (4)] with c0 or τsrc as the lone predictor (Figs. 3 and 4) show that we are able to reproduce ensemble member hindcast u¯ using M and xe. It remains to be shown that for a parameter value x′ not used in any ensemble member hindcast to fit and populate M, this M still accurately predicts the actual forecast model result using Eq. (4). To address this, we conduct single-parameter experiments for 1 February, 1 March, and 1 April 2014, generating a 30-day NAVGEM-HA hindcast ensemble for each initial date with only τsrc varying among members. For each case, a least squares solution to Eq. (4) produced an optimized τsrc value τ*. In each panel of Fig. 9, the blue curve is the u¯ profile estimated by evaluating Eq. (4) using τ* and the M for that initial date, while the green curve shows the corresponding u¯ profile from a new NAVGEM-HA hindcast run using this same GWD parameter setting. The agreement between the two is close in all three cases, even for the 1 April 2014 case (Fig. 9c) where τ* lies outside the range of values used to generate the original hindcast ensemble and to define M. The agreement is worst in the upper stratosphere, where the uncertainty arises in forecasting the stratospheric SAO rather than the QBO.

Fig. 9.

Profiles of equatorial u¯ from the NAVGEM-HA analysis (black), predicted optimal fit (blue), hindcast using optimal τsrc value (green), and hindcast using default τsrc value (red). Profiles are valid for 30-day lead times from initialization on (a) 1 Feb 2014, (b) 1 Mar 2014, and (c) 1 Apr 2014.

Fig. 9.

Profiles of equatorial u¯ from the NAVGEM-HA analysis (black), predicted optimal fit (blue), hindcast using optimal τsrc value (green), and hindcast using default τsrc value (red). Profiles are valid for 30-day lead times from initialization on (a) 1 Feb 2014, (b) 1 Mar 2014, and (c) 1 Apr 2014.

The red curves in Fig. 9 show corresponding NAVGEM-HA results using the default value of τsrc (4.5 mPa), and in February 2014 (Fig. 9a) the optimal τsrc value is coincidentally near this default value, so almost no additional skill is gained by optimizing τsrc, whereas in April (Fig. 9c) a large decrease in τsrc to ~2.2 mPa leads to a notable skill benefit relative to the default curve above about 7 hPa. In these initial tests, a least squares solution to Eq. (4) to obtain optimal parameters was obtained through a deep stratospheric layer (70–1 hPa). The results showed a tendency for the technique to select parameters that increase skill mostly in the upper stratosphere, where the default hindcast is less skillful and the ensemble spread is largest, which may end up reducing skill in hindcasts of the QBO in the lower stratosphere (e.g., Fig. 9c). In the next section, we will present results in which optimal parameters are derived via least squares solution to Eq. (4) for both the deep stratospheric layer (70–1 hPa) and a narrower pressure range focusing on the lower stratosphere (70–10 hPa).

To summarize, we have determined from the single- and multiparameter experiments that the best choice of predictive model M relating x′ and u¯ will use phase speed width c0 or source strength τsrc as a single predictor. We have validated the optimization methodology by demonstrating that our assumption that u¯ varies linearly with x′ is justifiable, that the linear model M can faithfully reproduce u¯' given x′ for the ensemble members, and that M can feasibly predict the behavior of the full forecast model at an extended lead time. We note that even when optimizing the GWD parameter subset {c0, τsrc}, we observe collinearity of the response profiles and an R2 value of 0.68 (not shown). Since we are unable to optimize this two-parameter subset simultaneously with sufficient accuracy, we therefore optimize each parameter individually. Both parameters have been shown to have connections to aspects of gravity wave dynamics within the real atmosphere, which make them appealing choices for optimization as physically relevant parameters. For example, Beres et al. (2004) showed that the phase speed width c0 of a spectrum of gravity waves generated by convection is linked to the horizontal and vertical structure of the convective source. It has also been suggested that the source momentum flux derived from τsrc could be linked quantitatively to aspects of the parameterized subgrid-scale convective activity, with resultant impacts on the QBO (Chun and Baik 2002; Beres et al. 2005; Kim et al. 2013; Schirber et al. 2014; Kang et al. 2018). Equatorial wave forcing has also been shown to peak near equinox (Kim and Chun 2015), meaning that such convective GWD schemes could additionally serve to supplement underrepresented resolved wave forcing in a model.

4. Results

Accelerations due to parameterized wave–mean flow interaction usually peak below QBO shear zones as varying background zonal winds present critical levels, below which waves saturate and deposit momentum (e.g., McCormack et al. 2015; Anstey et al. 2016). Thus, the shape of the response profiles in M will be contingent on the shape of the initial u¯ profile, implying that M will vary with latitude, altitude range, and season. M also has an implicit dependence on lead time Δt. From Eq. (4), it follows that the optimized GWD source parameter x′ will share these dependencies, therefore we choose to analyze ensembles of hindcasts following the approach detailed previously for multiple latitudes, column depths, initial times, and lead times.

Results presented in this section will cover the NAVGEM-HA analysis years 2014, 2015, and 2016, sampled monthly by launching a hindcast ensemble at 0000 UTC on the first day of every month, for both c0 and τsrc single-predictor models. Each parameter will be optimized for 30-day and 100-day lead times and for column depths of 70–10 and 70–1 hPa. In addition to optimizing the equatorial source parameters, we will also explore the effect of varying τsrc in the Northern Hemisphere (NH) subtropics (20°N) and extratropics (60°N) on forecast skill both at that latitude and at the equator. We reiterate that these results may in part leverage the nonorographic GWD parameterization to compensate for model error due to, for example, insufficient representation of resolved wave forcing.

a. Subseasonal (30 day) optimization

Time series of optimal source strength τsrc values that minimize 30-day lead time AF error are summarized in Fig. 10b in comparison to analyzed equatorial stratospheric winds (Fig. 10a). The most striking feature is a semiannual variation in optimal τsrc for the 70–10-hPa column depth with equinoctial peaks. Spectral analysis indicates that the semiannual frequency is indeed the strongest spectral mode (not shown). This result is consistent with the hypothesis that source strength is related to underlying convective activity, which peaks at equinox in the equatorial region. Since large-scale equatorial waves are also generated by convection, the optimization could simultaneously be leveraging the nonorographic GWD to correct a resolved momentum budget deficit. During the anomalous so-called QBO disruption period of 2016 (Osprey et al. 2016), the equinoctial peaks vanish. Optimal τsrc values at this time may not be trustworthy considering both the difficulty in hindcasting the event and that the anomalous drag in this event came from an extratropical source, rather than from tropical GWD (Osprey et al. 2016; Coy et al. 2017; Barton and McCormack 2017).

Fig. 10.

Optimization of source strength τsrc and phase speed width c0. (a) NAVGEM-HA analyzed zonal-mean zonal wind over the equator (2°S–2°N average). (b) Values of τsrc for the default setting (black) and minimizing 30-day hindcast error over depths of 70–10 hPa (green) and 70–1 hPa (orange) as a function of initialization date. (c) RMS error of optimal hindcast (solid) and R2 measure of M (dotted) with color corresponding to vertical depth as in (b). Lines dashed with black are RMS error of the default hindcast calculated over the depth corresponding to the secondary color. (d), (e) As in (b) and (c), but for phase speed width c0.

Fig. 10.

Optimization of source strength τsrc and phase speed width c0. (a) NAVGEM-HA analyzed zonal-mean zonal wind over the equator (2°S–2°N average). (b) Values of τsrc for the default setting (black) and minimizing 30-day hindcast error over depths of 70–10 hPa (green) and 70–1 hPa (orange) as a function of initialization date. (c) RMS error of optimal hindcast (solid) and R2 measure of M (dotted) with color corresponding to vertical depth as in (b). Lines dashed with black are RMS error of the default hindcast calculated over the depth corresponding to the secondary color. (d), (e) As in (b) and (c), but for phase speed width c0.

Unlike the lower stratosphere fit, the whole stratosphere fit (70–1 hPa, orange curve) has a much more pronounced annual cycle, with a large primary peak in boreal winter and a weaker secondary maximum near July. As discussed previously, large AF error in the upper stratosphere would overwhelm the relatively skillful QBO prediction in the lower stratosphere, yielding results incorrectly synchronized with the easterly phase of the SAO. However, because the hindcast error is not due primarily to insufficient gravity wave driving, the best fit RMS error for the whole stratosphere fit is on average 5.3 m s−1 higher than that of the lower stratosphere fit (see Fig. 10c). We note that the R2 measure of model fidelity is over 0.9 for most of the analysis period in both column depth cases.

The optimal values of phase speed width in Fig. 10d reveal a less coherent time evolution than for τsrc. The annual cycle is more significant for the 70–10-hPa fit. More analysis years are required to discern a characteristic pattern in optimal c0 as a function of initial condition for 30-day forecasts. In the whole stratosphere case, we find that the time evolution behaves similarly to that of τsrc, which again is a by-product of interaction with SAO easterly phases in the tropical upper stratosphere. Due to larger ensemble spread in the upper stratosphere (cf. Figs. 3 and 4), we expected that the best fit c0 hindcasts would have lower RMS error than their τsrc counterparts, which we find to be the case. We note that the inability to hindcast the easterly SAO phase is more conspicuous in Fig. 10e, where there are significant peaks in RMS error during solstice.

b. Seasonal (100 day) optimization

In Fig. 11b, we present similar τsrc time series as Fig. 10b, calculated here for a lead time of 100 days. While the optimal τsrc value varied semiannually for 30-day hindcasts, the green curve in Fig. 11b indicates ostensibly that τsrc varies with QBO phase, with higher values during the westerly phase. In the standard QBO paradigm, the westerly phase is thought to be maintained primarily by resolved waves (Giorgetta et al. 2006). Thus, we hypothesize that more parameterized gravity wave forcing is required to counterbalance the tendency of the model winds to trend easterly over longer forecast periods. Another possibility is that hindcasts of a QBO transition period require more GWD contribution to lower the QBO shear zone. In this case, our assumption of a τsrc value that does not vary during the hindcast may be especially inaccurate when optimizing long hindcasts that include a transition period. More analysis years are needed to substantiate these hypotheses. There is some weak semiannual activity suggested by the curve, but we note that the 100-day hindcast length is near the resolvable Nyquist limit of the 180-day semiannual period (90 days), which would cause the semiannual variability in the optimal τsrc curve to be artificially weakened. The whole stratosphere fit (orange curve) has no remarkable characteristics. As expected, the fidelity of the model M as measured by R2 decreased as the lead time increased (see Fig. 11c). However, the lower stratosphere fit (green curve) continues to exhibit a relatively high R2 value over the analysis period. Difficulty in hindcasting the QBO disruption at seasonal lead time can be easily seen in the increasing RMS error values in early 2016.

Fig. 11.

As in Fig. 10, but for minimization of 100-day hindcast error.

Fig. 11.

As in Fig. 10, but for minimization of 100-day hindcast error.

The lack of a similar peak in RMS error in the 30-day optimization is indicative of the relatively long radiative relaxation time scales in the lower stratosphere, where a forecast of the QBO can be nominally “good” for some time simply by virtue of initializing the model with accurate QBO initial conditions, which can then persist in the forecasts (but not descend realistically) for many days absent other forcings. This is also reflected in the contrast between green curves in Figs. 10b and 11b, where we are attempting to optimize forecasts in different time regimes. In the case of subseasonal (i.e., 30 day) prediction, in which we expect there is some intrinsic skill in forecasting QBO winds, we find that the convectively generated wave spectrum must be supported by a parameterized component (i.e., the nonorographic GWD scheme), which is realized by a semiannually varying source function. On the other hand, errors in prediction at the seasonal scale are due more to inadequate simulation of the QBO, so the gravity wave drag is optimized to compensate, which leads to a source function synchronized with the QBO. It is possible that longer free-running simulations of the QBO would necessitate a source function that varies both semiannually and with the QBO. In general, we have found that a constant GWD source strength in the equatorial band is not optimal for minimizing hindcast errors in QBO winds.

Results for the phase speed width optimization for a 100-day lead time (Fig. 11d) are similar in many ways to those of τsrc in Fig. 11b, particularly in the temporal evolution of optimal c0 curves. The best fit RMS error is as low as in the subseasonal case (2–4 m s−1) during the easterly phase of the QBO (latter half of 2014), a feature not found in the τsrc optimization. In terms of RMS error, c0 is more effective at increasing forecast skill than τsrc at both subseasonal and seasonal time scales; however, it also suffers from lower R2.

c. Subtropical and extratropical optimization

We have also explored the use of our optimization technique at higher latitudes. The main caveat is that, as shown in Fig. 12, the long-range predictive skill of our linear model Eq. (4) in optimizing GWD for zonal-mean zonal winds degrades with latitude, as quantified by calculating the R2 measure as a function of lead time for ensembles run at various latitudes. The results for the equator, 20°, and 60°N are shown in Fig. 12, where we have again minimized the AF error at these latitudes for model levels from 70 to 10 hPa. We see that for lead times from 10 to 100 days, GWD parameters can be accurately optimized to reproduce analyzed zonal-mean zonal winds over the equator (Fig. 12a). Boreal summer months are the most difficult to model linearly in the tropics. On the other hand, predictive skill at 20° and 60°N degrades quickly as forecast lead time increases. For 10-day forecasts, the subtropics (Fig. 12b) hold some amount of skill (R2 > 0.9), but this value is not sustained for longer lead times. Boreal summer months are the easiest to model linearly, in contrast to the result at the equator. In the extratropics (Fig. 12c), only a handful of months show reasonable skill at 10 days, and there are very few cases of R2 > 0.6 for any month at any longer lead time.

Fig. 12.

The R2 measures of model hindcast fidelity as a function of lead time for (a) the equator, (b) 20°N, and (c) 60°N, averaged by month (color coded). The R2 is calculated for ensembles with varying τsrc launched from the first of every month for NAVGEM-HA analysis years 2014, 2015, and 2016, as in previous sections.

Fig. 12.

The R2 measures of model hindcast fidelity as a function of lead time for (a) the equator, (b) 20°N, and (c) 60°N, averaged by month (color coded). The R2 is calculated for ensembles with varying τsrc launched from the first of every month for NAVGEM-HA analysis years 2014, 2015, and 2016, as in previous sections.

It is perhaps unsurprising that the methodology as formulated for the tropical stratosphere is ineffective at minimizing error in the extratropics at long lead times. In the tropics, wind and temperature tendencies are driven locally by GWD, while in the extratropics, the majority of nonorographic wave breaking occurs in the mesosphere and above, and influence on the stratosphere is exerted through nonlocal effects such as downward control. The extratropical stratospheric circulation is also driven more directly by planetary wave drag and orographic GWD, especially in local winter (Alexander and Rosenlof 2003). As such, the influence of our parameterized nonorographic GWD in this region may be small, and so we cannot expect to effectively minimize the hindcast error without building downward control effects of upper-level GWD and influences of orographic GWD into our optimization technique. In short, a careful redesign of the cost function is needed to optimize the GWD in the extratropics.

Although we have found that optimizing GWD at higher latitudes by fitting against analysis wind at that latitude is not effective for lead times longer than about 10 days, we additionally attempted to fit against analysis wind at the equator to observe effects on the simulated QBO due to changes in GWD at higher latitudes. Specifically, we were interested in gaining hindcast skill for the QBO disruption event by altering the sub- and extratropical waveguide through the parameterized GWD. Barton and McCormack (2017) suggested that a serendipitous reduction in easterlies in the Northern Hemisphere subtropics allowed the anomalous extratropical flux to reach the equator and induce the disruption. Thus, we explore whether changes in sub- and extratropical winds as a result of changes in the GWD source strength can affect hindcast winds at the equator.

We find that there is very little effect on the equatorial hindcast wind profile in the lower stratosphere by changing τsrc at 20° or 60°N. Figure 13 shows ensemble hindcast u¯ profiles for 30-day hindcasts launched 1 February 2016 and 100-day hindcasts launched 1 December 2015. In both hindcast sets, the final valid date is in early March when the easterly disruption had emerged in the analysis. For 30-day hindcasts, changing τsrc at 20° or 60°N has almost no observable effect on the winds below about 7 hPa, and the ensemble spread is only slightly larger for the 100-day hindcasts. In the upper stratosphere we are able to see some variability in equatorial winds, but in each plot except Fig. 13a the winds vary nonlinearly with τsrc. In no case does the ensemble spread appear to be able to account for the wedge of easterlies that developed near 40 hPa. Lin et al. (2019) found that the anomalous forcing leading to the event was not in the zonal mean but localized in longitude by single wave packet. It may therefore be more suitable for optimization of GWD source parameters during the QBO disruption period if our cost function was redefined to fit the zonal wind in the 2D longitude–height plane, rather than fitting the zonal-mean zonal wind in height alone.

Fig. 13.

Profiles of equatorial u¯ valid for early March 2016, for (top) 30- and (bottom) 100-day ensembles launched with τsrc varying at (a),(d) the equator, (b),(e) 20°N, and (c),(f) 60°N. The analysis profile valid for the final hindcast date is shown in black, and the ensemble member profiles are color-coded by τsrc value. For the 20 and 60°N optimizations, the equatorial τsrc was given its default value of 4.5 mPa.

Fig. 13.

Profiles of equatorial u¯ valid for early March 2016, for (top) 30- and (bottom) 100-day ensembles launched with τsrc varying at (a),(d) the equator, (b),(e) 20°N, and (c),(f) 60°N. The analysis profile valid for the final hindcast date is shown in black, and the ensemble member profiles are color-coded by τsrc value. For the 20 and 60°N optimizations, the equatorial τsrc was given its default value of 4.5 mPa.

5. Implications for prediction

The motivation for studying years of analysis with this technique is to inform decisions about how the parameterized GWD source in a particular model should be prescribed in order to maximize predictive skill. As we are unable to verify forecasts of future events, we require a robust characterization of the optimal source parameters so that we can be confident our forecasts will be accurate for a given season. The challenge lies in distilling the information yielded by the least squares optimization solutions to the linearized version of the problem in Eq. (4) (e.g., Fig. 10b) into a meaningful characterization of optimal τsrc for true forecasts rather than hindcasts. A physical link between the atmospheric state and the source parameters is a desirable feature that will ensure the nonorographic GWD scheme can dynamically adjust to the optimal settings as the forecast progresses. For example, the implied relationship between optimal τsrc and convection evident in the semiannual variability in optimized GWD parameters (e.g., Fig. 10b) validates the future use of a convection-based GWD parameterization, and provides constraints and guidelines for how the convection parameterization should modulate GW momentum and phase speeds at tropical source levels in the model.

Following the results of section 4a, our GWD parameterization is optimized by introducing an explicit semiannual variability to the source strength function at the equator that better reproduces zonal-mean zonal winds in the tropical lower stratosphere. Figure 14 illustrates this result by comparing the observed QBO behavior (Fig. 14a) with the time series of τsrc values imposed in the NAVGEM-HA GWD parameterization in various 30-day hindcast experiments over the January 2014–December 2016 period (Fig. 14b) and the corresponding hindcast RMS error (Fig. 14c). The semiannual component of the green curve in Fig. 10b, which we determined to be the primary mode of variability through spectral analysis, is extracted and plotted in purple against the original in Fig. 14b. Indeed, this single mode accounts for a significant fraction of variance in the full time series, except during the QBO disruption event in early 2016. Due to its simplicity, a regular semiannual variation in τsrc is easily implemented in our scheme using the amplitude and phase values obtained from the spectral analysis. We redefine the equatorial source strength as a semiannual (period = one-half year) sine function of day with mean 2.285 mPa, amplitude 1.191 mPa, and phase −1.86 radians. This means that τsrc will vary over the course of the hindcast, unlike in the previous experiments where τsrc was constant. The yearly peaks of this function occur in early April and October, immediately following equinox. For comparison, we also show the results of keeping τsrc constant but lowering it to the mean value of the optimal curve (2.285 mPa, orange curve).

Fig. 14.

(a) NAVGEM-HA analysis zonal-mean zonal wind over the equator (2°S–2°N average). (b) Comparison of τsrc values: default (black), optimized for 30-day hindcast over 70–10 hPa (green, as in Fig. 10b), mean of the optimized values (orange), and the semiannual component with mean of the optimized values (purple). (c) RMS error of 30-day hindcasts with τsrc values as in (b).

Fig. 14.

(a) NAVGEM-HA analysis zonal-mean zonal wind over the equator (2°S–2°N average). (b) Comparison of τsrc values: default (black), optimized for 30-day hindcast over 70–10 hPa (green, as in Fig. 10b), mean of the optimized values (orange), and the semiannual component with mean of the optimized values (purple). (c) RMS error of 30-day hindcasts with τsrc values as in (b).

The RMS errors of hindcasts using default, best fit, best fit mean, and semiannually varying τsrc are compared in Fig. 14c. Lowering the value of τsrc to the mean of the best fit curve reduces RMS error in most cases. By introducing semiannual variability to the equatorial source strength function, we further reduce RMS error by an average of 1.5 m s−1 compared to using the default τsrc value of 4.5 mPa. The difference between the best fit mean (orange) and semiannual τsrc (purple) curves is statistically significant at the 90% confidence level. Additionally, we find only marginal difference in using a simple semiannual function of τsrc compared to the actual optimal τsrc values, which gives us confidence that forecast errors of QBO winds will be reduced also for predictions made outside of the 2014–16 analysis period shown, given that this variability is a robust feature of the source strength optimized for subseasonal prediction in any typical QBO season. The QBO disruption period shows poor predictability using the semiannual τsrc because the semiannual function diverges significantly from the optimal value.

Performing these subseasonal and seasonal optimizations within a longer analysis period spanning many QBO cycles may provide evidence that other frequency modes are also significant in the time evolution of optimal GWD source parameter values. For example, the largest values of τsrc in the Fig. 14b green curve occur during a transition period of the QBO, which corresponds with our hypothesis that the QBO period was the dominant mode of variability in seasonal τsrc optimization (Fig. 11b). This is evidence that long-term effects of the parameterized GWD (its role in the descent of QBO shear zones, for example) have an influence on shorter forecasts. These modes of variability could also be easily introduced but could conceivably require some awareness of the forecast length if τsrc is kept constant over the forecast length rather than allowed to vary. Improvements to the resolved wave spectrum through updated convection parameterization or model resolution, for example, may further shift the role of the nonorographic GWD scheme to a more passive one as the model becomes less reliant on the parameterization to provide missing momentum.

6. Conclusions

We have formulated a methodology for determining the values of GWD source parameters that minimize AF error within a domain of interest for a desired forecast lead time. The utilization of this technique is straightforward, requiring only ensembles of hindcasts with source parameters varying among ensemble members. Least squares minimization reveals the effect of perturbing a source parameter on the final hindcast zonal winds in the target domain, from which we can calculate the parameter value that brings the hindcast wind closest to the analysis state. This simultaneously corrects for parameterized GWD error as well as error due to insufficient resolved wave forcing without explicitly calculating and compensating for a momentum budget deficit. Applying the technique, we have found that optimizing a single source parameter works best, and that source momentum flux τsrc and phase speed spectrum width c0 are the best choices as predictors. Both are effective in reproducing ensemble member hindcast profiles of QBO winds over the equator using a simple linear model, even at extended lead time.

Analysis of the 2014–16 period suggests that values of τsrc and c0 that minimize forecast error in the equatorial lower stratosphere are a function of initial condition, rather than static values as they are currently in NAVGEM-HA. Interestingly, the characteristics of the evolution vary between subseasonal and seasonal forecast lengths. For 30-day forecasts, there is a prominent semiannual variability in τsrc and an apparent secondary signal synchronized with the QBO, while only this secondary signal manifests in 100-day forecasts. The GWD parameterization appears to become more important in time as the forecast loses memory of the initial state, which is more relevant for QBO simulation in longer, seasonal forecasts than in shorter forecasts. In our attempt to optimize subseasonal prediction, we have drawn a connection between τsrc and convective activity, which is substantiated by previous work deriving source momentum flux from convection (e.g., Beres et al. 2004; Richter et al. 2010; Schirber et al. 2014). The time series of optimal c0 parallels that of τsrc, except that a significant annual cycle also exists for the 30-day lead time.

We reiterate that these results are specific to this configuration of NAVGEM-HA and will likely differ with results from applying the procedure to other models. For example, models with higher resolution or more skillful forecasts of primary QBO drivers may find that the necessary forcing supplied by the nonorographic GWD parameterization is less than in our results.

A predictability time scale of 10 days or less in the sub- and extratropical Northern Hemisphere, combined with additional complications of downward control and parameterized orographic GWD, precludes us from identifying optimal source parameter values at those latitudes for subseasonal or seasonal prediction. We also found that for both 30- and 100-day hindcasts, there was no effect on the hindcast QBO winds at the equator by changing the GWD source at higher latitudes. Thus, we were unable to improve hindcasts of the 2016 QBO disruption through GWD parameterization with optimal source momentum flux.

However, we found that replacing the constant equatorial τsrc value in NAVGEM-HA with one that varies semiannually with season resulted in reduced RMS error in stratospheric equatorial zonal-mean zonal wind throughout the analysis period. This simple implementation was informed by our analysis of optimal τsrc (see Fig. 10) and the satisfactory result is evidence that a more dynamic parameterization of GW source-level parameters in the equatorial tropics is necessary to capture a realistic source of the GW momentum that is subsequently deposited in the lower stratosphere and drives a QBO in the model. Because the semiannual mode accounts for the majority of variance in the time series, we are confident that the change to equatorial τsrc will be beneficial for true forecasts as well as hindcast experiments.

With a priori knowledge that τsrc could vary semiannually, it is also possible for us to impose the semiannual variability in the model and optimize for the amplitude or mean of the oscillation. This has implications especially for optimization of longer hindcasts, where the difference between an assumed constant τsrc value (as in section 4) and a value that varies with day of year may be significant. Future studies will focus on developing a more detailed and accurate function for τsrc that could be determined by optimizing decades of analysis winds encompassing many QBO cycles.

Acknowledgments

The authors acknowledge support from the Chief of Naval Research. This research was performed while Cory Barton held an NRC Research Associateship award at the Naval Research Laboratory. NAVGEM-HA analyses and hindcasts were produced under a grant of computer time from the Department of Defense High Performance Computing Modernization Program, and the fields used in this study are available upon request by contacting the authors.

REFERENCES

REFERENCES
Alexander
,
M. J.
, and
J. R.
Holton
,
1997
:
A model study of zonal forcing in the equatorial stratosphere by convectively induced gravity waves
.
J. Atmos. Sci.
,
54
,
408
419
, https://doi.org/10.1175/1520-0469(1997)054<0408:AMSOZF>2.0.CO;2.
Alexander
,
M. J.
, and
K. H.
Rosenlof
,
2003
:
Gravity-wave forcing in the stratosphere: Observational constraints from the Upper Atmosphere Research Satellite and implications for parameterization in global models
.
J. Geophys. Res.
,
108
,
4597
, https://doi.org/10.1029/2003JD003373.
Alexander
,
M. J.
, and Coauthors
,
2010
:
Recent developments in gravity-wave effects in climate models and the global distribution of gravity-wave momentum flux from observations and models
.
Quart. J. Roy. Meteor. Soc.
,
136
,
1103
1124
, https://doi.org/10.1002/QJ.637.
Allen
,
S. J.
, and
R. A.
Vincent
,
1995
:
Gravity wave activity in the lower atmosphere: Seasonal and latitudinal variations
.
J. Geophys. Res.
,
100
,
1327
1350
, https://doi.org/10.1029/94JD02688.
Anstey
,
J. A.
,
J. F.
Scinocca
, and
M.
Keller
,
2016
:
Simulating the QBO in an atmospheric general circulation model: Sensitivity to resolved and parameterized forcing
.
J. Atmos. Sci.
,
73
,
1649
1665
, https://doi.org/10.1175/JAS-D-15-0099.1.
Baldwin
,
M. P.
, and Coauthors
,
2001
:
The quasi-biennial oscillation
.
Rev. Geophys.
,
39
,
179
229
, https://doi.org/10.1029/1999RG000073.
Barton
,
C. A.
, and
J. P.
McCormack
,
2017
:
Origin of the 2016 QBO disruption and its relationship to extreme El Niño events
.
Geophys. Res. Lett.
,
44
,
11 150
11 157
, https://doi.org/10.1002/2017GL075576.
Beres
,
J. H.
,
M. J.
Alexander
, and
J. R.
Holton
,
2004
:
A method of specifying the gravity wave spectrum above convection based on latent heating properties and background wind
.
J. Atmos. Sci.
,
61
,
324
337
, https://doi.org/10.1175/1520-0469(2004)061<0324:AMOSTG>2.0.CO;2.
Beres
,
J. H.
,
R. R.
Garcia
,
B. A.
Boville
, and
F.
Sassi
,
2005
:
Implementation of a gravity wave source spectrum parameterization dependent on the properties of convection in the Whole Atmosphere Community Climate Model (WACCM)
.
J. Geophys. Res.
,
110
,
D10108
, https://doi.org/10.1029/2004JD005504.
Bushell
,
A. C.
,
N.
Butchart
,
S. H.
Derbyshire
,
D. R.
Jackson
,
G. J.
Shutts
,
S. B.
Vosper
, and
S.
Webster
,
2015
:
Parameterized gravity wave momentum fluxes from sources related to convection and large-scale precipitation processes in a global atmosphere model
.
J. Atmos. Sci.
,
72
,
4349
4371
, https://doi.org/10.1175/JAS-D-15-0022.1.
Butchart
,
N.
, and Coauthors
,
2018
:
Overview of experiment design and comparison of models participating in phase 1 of the SPARC Quasi-Biennial Oscillation initiative (QBOi)
.
Geosci. Model Dev.
,
11
,
1009
1032
, https://doi.org/10.5194/gmd-11-1009-2018.
Chen
,
X.
,
O. M.
Pauluis
, and
F.
Zhang
,
2018
:
Regional simulation of Indian summer monsoon intraseasonal oscillations at gray-zone resolution
.
Atmos. Chem. Phys.
,
18
,
1003
1022
, https://doi.org/10.5194/acp-18-1003-2018.
Chun
,
H.-Y.
, and
J.-J.
Baik
,
2002
:
An updated parameterization of convectively forced gravity wave drag for use in large-scale models
.
J. Atmos. Sci.
,
59
,
1006
1017
, https://doi.org/10.1175/1520-0469(2002)059<1006:AUPOCF>2.0.CO;2.
Coy
,
L.
,
P. A.
Newman
,
S.
Pawson
, and
L. R.
Lait
,
2017
:
Dynamics of the disrupted 2015/16 quasi-biennial oscillation
.
J. Climate
,
30
,
5661
5674
, https://doi.org/10.1175/JCLI-D-16-0663.1.
de la Cámara
,
A.
,
F.
Lott
, and
M.
Abalos
,
2016
:
Climatology of the middle atmosphere in LMDz: Impact of source-related parameterizations of gravity wave drag
.
J. Adv. Model. Earth Syst.
,
8
,
1507
1525
, https://doi.org/10.1002/2016MS000753.
Dunkerton
,
T. J.
,
1982
:
Theory of the mesopause semiannual oscillation
.
J. Atmos. Sci.
,
39
,
2681
2690
, https://doi.org/10.1175/1520-0469(1982)039<2681:TOTMSO>2.0.CO;2.
Dunkerton
,
T. J.
,
1997
:
The role of gravity waves in the quasi-biennial oscillation
.
J. Geophys. Res.
,
102
,
26 053
26 076
, https://doi.org/10.1029/96JD02999.
Eckermann
,
S. D.
,
2011
:
Explicitly stochastic parameterization of nonorographic gravity wave drag
.
J. Atmos. Sci.
,
68
,
1749
1765
, https://doi.org/10.1175/2011JAS3684.1.
Eckermann
,
S. D.
, and Coauthors
,
2009
:
High-altitude data assimilation system experiments for the northern summer mesosphere season of 2007
.
J. Atmos. Sol.-Terr. Phys.
,
71
,
531
551
, https://doi.org/10.1016/j.jastp.2008.09.036.
Eckermann
,
S. D.
, and Coauthors
,
2018
:
High-altitude (0–100 km) global atmospheric reanalysis system: Description and application to the 2014 austral winter of the Deep Propagating Gravity Wave Experiment (DEEPWAVE)
.
Mon. Wea. Rev.
,
146
,
2639
2666
, https://doi.org/10.1175/MWR-D-17-0386.1.
Fritts
,
D. C.
, and
M. J.
Alexander
,
2003
:
Gravity wave dynamics and effects in the middle atmosphere
.
Rev. Geophys.
,
41
,
1003
, https://doi.org/10.1029/2001RG000106.
Garcia
,
R. R.
, and
B. A.
Boville
,
1994
:
“Downward control” of the mean meridional circulation and temperature distribution of the polar winter stratosphere
.
J. Atmos. Sci.
,
51
,
2238
2245
, https://doi.org/10.1175/1520-0469(1994)051<2238:COTMMC>2.0.CO;2.
Garcia
,
R. R.
,
T. J.
Dunkerton
,
R. S.
Lieberman
, and
R.
Vincent
,
1997
:
Climatology of the semiannual oscillation of the tropical middle atmosphere
.
J. Geophys. Res.
,
102
,
26 019
26 032
, https://doi.org/10.1029/97JD00207.
Garcia
,
R. R.
,
D. R.
Marsh
,
D. E.
Kinnison
,
B. A.
Boville
, and
F.
Sassi
,
2007
:
Simulation of secular trends in the middle atmosphere, 1950-2003
.
J. Geophys. Res.
,
112
,
D09301
, https://doi.org/10.1029/2006JD007485.
Garcia
,
R. R.
,
A. K.
Smith
,
D. E.
Kinnison
,
Á.
Cámara
, and
D. J.
Murphy
,
2017
:
Modification of the gravity wave parameterization in the Whole Atmosphere Community Climate Model: Motivation and results
.
J. Atmos. Sci.
,
74
,
275
291
, https://doi.org/10.1175/JAS-D-16-0104.1.
Garfinkel
,
C. I.
, and
D. L.
Hartmann
,
2011
:
The influence of the quasi-biennial oscillation on the troposphere in winter in a hierarchy of models. Part I: Simplified dry GCMs
.
J. Atmos. Sci.
,
68
,
1273
1289
, https://doi.org/10.1175/2011JAS3665.1.
Garfinkel
,
C. I.
,
C.
Schwartz
,
D. I.
Domeisen
,
S.-W.
Son
,
A. H.
Butler
, and
I. P.
White
,
2018
:
Extratropical atmospheric predictability from the quasi-biennial oscillation in subseasonal forecast models
.
J. Geophys. Res. Atmos.
,
123
,
7855
7866
, https://doi.org/10.1029/2018JD028724.
Geller
,
M. A.
, and Coauthors
,
2016
:
Modeling the QBO—Improvements resulting from higher-model vertical resolution
.
J. Adv. Model. Earth Syst.
,
8
,
1092
1105
, https://doi.org/10.1002/2016MS000699.
Giorgetta
,
M. A.
,
E.
Manzini
, and
E.
Roeckner
,
2002
:
Forcing of the quasi-biennial oscillation from a broad spectrum of atmospheric waves
.
Geophys. Res. Lett.
,
29
,
1245
, https://doi.org/10.1029/2002GL014756.
Giorgetta
,
M. A.
,
E.
Manzini
,
E.
Roeckner
,
M.
Esch
, and
L.
Bengtsson
,
2006
:
Climatology and forcing of the quasi-biennial oscillation in the MAECHAM5 model
.
J. Climate
,
19
,
3882
3901
, https://doi.org/10.1175/JCLI3830.1.
Hamilton
,
K.
,
R. J.
Wilson
, and
R. S.
Hemler
,
1999
:
Middle atmosphere simulated with high vertical and horizontal resolution versions of a GCM: Improvements in the cold pole bias and generation of a QBO-like oscillation in the tropics
.
J. Atmos. Sci.
,
56
,
3829
3846
, https://doi.org/10.1175/1520-0469(1999)056<3829:MASWHV>2.0.CO;2.
Haynes
,
P. H.
,
M. E.
McIntyre
,
T. G.
Shepherd
,
C. J.
Marks
, and
K. P.
Shine
,
1991
:
On the “downward control” of extratropical diabatic circulations by eddy-induced mean zonal forces
.
J. Atmos. Sci.
,
48
,
651
678
, https://doi.org/10.1175/1520-0469(1991)048<0651:OTCOED>2.0.CO;2.
Hitchman
,
M. H.
, and
A. S.
Huesmann
,
2009
:
Seasonal influence of the quasi-biennial oscillation on stratospheric jets and Rossby wave breaking
.
J. Atmos. Sci.
,
66
,
935
946
, https://doi.org/10.1175/2008JAS2631.1.
Hitchman
,
M. H.
,
J. C.
Gille
,
C. D.
Rodgers
, and
G.
Brasseur
,
1989
:
The separated polar winter stratopause: A gravity wave driven climatological feature
.
J. Atmos. Sci.
,
46
,
410
422
, https://doi.org/10.1175/1520-0469(1989)046<0410:TSPWSA>2.0.CO;2.
Hogan
,
T.
, and Coauthors
,
2014
:
The Navy Global Environmental Model
.
Oceanography
,
27
(
3
),
116
125
, https://doi.org/10.5670/oceanog.2014.73.
Holt
,
L. A.
,
M. J.
Alexander
,
L.
Coy
,
A.
Molod
,
W.
Putman
, and
S.
Pawson
,
2016
:
Tropical waves and the quasi-biennial oscillation in a 7-km global climate simulation
.
J. Atmos. Sci.
,
73
,
3771
3783
, https://doi.org/10.1175/JAS-D-15-0350.1.
Holton
,
J. R.
,
1983
:
The influence of gravity wave breaking on the general circulation of the middle atmosphere
.
J. Atmos. Sci.
,
40
,
2497
2507
, https://doi.org/10.1175/1520-0469(1983)040<2497:TIOGWB>2.0.CO;2.
Holton
,
J. R.
, and
R. S.
Lindzen
,
1972
:
An updated theory for the quasi-biennial cycle of the tropical stratosphere
.
J. Atmos. Sci.
,
29
,
1076
1080
, https://doi.org/10.1175/1520-0469(1972)029<1076:AUTFTQ>2.0.CO;2.
Holton
,
J. R.
, and
H.-C.
Tan
,
1980
:
The influence of the equatorial quasi-biennial oscillation on the global circulation at 50 mb
.
J. Atmos. Sci.
,
37
,
2200
2208
, https://doi.org/10.1175/1520-0469(1980)037<2200:TIOTEQ>2.0.CO;2.
Horinouchi
,
T.
, and Coauthors
,
2003
:
Tropical cumulus convection and upward-propagating waves in middle-atmosphere GCMs
.
J. Atmos. Sci.
,
60
,
2765
2782
, https://doi.org/10.1175/1520-0469(2003)060<2765:TCCAUW>2.0.CO;2.
Janiga
,
M. A.
,
C. J.
Schreck
,
J. A.
Ridout
,
M.
Flatau
,
N. P.
Barton
,
E. J.
Metzger
, and
C. A.
Reynolds
,
2018
:
Subseasonal forecasts of convectively coupled equatorial waves and the MJO: Activity and predictive skill
.
Mon. Wea. Rev.
,
146
,
2337
2360
, https://doi.org/10.1175/MWR-D-17-0261.1.
Kang
,
M.-J.
,
H.-Y.
Chun
,
Y.-H.
Kim
,
P.
Preusse
, and
M.
Ern
,
2018
:
Momentum flux of convective gravity waves derived from an offline gravity wave parameterization. Part II: Impacts on the quasi-biennial oscillation
.
J. Atmos. Sci.
,
75
,
3753
3775
, https://doi.org/10.1175/JAS-D-18-0094.1.
Kawatani
,
Y.
,
K.
Sato
,
T. J.
Dunkerton
,
S.
Watanabe
,
S.
Miyahara
, and
M.
Takahashi
,
2010
:
The roles of equatorial trapped waves and internal inertia–gravity waves in driving the quasi-biennial oscillation. Part I: Zonal mean wave forcing
.
J. Atmos. Sci.
,
67
,
963
980
, https://doi.org/10.1175/2009JAS3222.1.
Kim
,
Y.-H.
, and
H.-Y.
Chun
,
2015
:
Contributions of equatorial wave modes and parameterized gravity waves to the tropical QBO in HadGEM2
.
J. Geophys. Res. Atmos.
,
120
,
1065
1090
, https://doi.org/10.1002/2014JD022174.
Kim
,
Y.-H.
,
A. C.
Bushell
,
D. R.
Jackson
, and
H.-Y.
Chun
,
2013
:
Impacts of introducing a convective gravity-wave parameterization upon the QBO in the Met Office United Model
.
Geophys. Res. Lett.
,
40
,
1873
1877
, https://doi.org/10.1002/grl.50353.
Kim
,
Y.-J.
,
S. D.
Eckermann
, and
H.-Y.
Chun
,
2003
:
An overview of the past, present, and future of gravity-wave drag parameterization for numerical climate and weather prediction models
.
Atmos.–Ocean
,
41
,
65
98
, https://doi.org/10.3137/ao.410105.
Krismer
,
T. R.
,
M. A.
Giorgetta
,
J. S.
von Storch
, and
I.
Fast
,
2015
:
The influence of the spectral truncation on the simulation of waves in the tropical stratosphere
.
J. Atmos. Sci.
,
72
,
3819
3828
, https://doi.org/10.1175/JAS-D-14-0240.1.
Kuhl
,
D. D.
,
T. E.
Rosmond
,
C. H.
Bishop
,
J.
McLay
, and
N. L.
Baker
,
2013
:
Comparison of hybrid ensemble/4DVar and 4DVar within the NAVDAS-AR data assimilation framework
.
Mon. Wea. Rev.
,
141
,
2740
2758
, https://doi.org/10.1175/MWR-D-12-00182.1.
Lin
,
P.
,
I.
Held
, and
Y.
Ming
,
2019
:
The early development of the 2015/2016 quasi-biennial oscillation disruption
.
J. Atmos. Sci.
,
76
,
821
836
, https://doi.org/10.1175/JAS-D-18-0292.1.
Lindzen
,
R. S.
,
1981
:
Turbulence and stress owing to gravity wave and tidal breakdown
.
J. Geophys. Res.
,
86
,
9707
9714
, https://doi.org/10.1029/JC086iC10p09707.
Lindzen
,
R. S.
, and
J. R.
Holton
,
1968
:
A theory of the quasi-biennial oscillation
.
J. Atmos. Sci.
,
25
,
1095
1107
, https://doi.org/10.1175/1520-0469(1968)025<1095:ATOTQB>2.0.CO;2.
Long
,
D. J.
,
D. R.
Jackson
, and
J.
Thuburn
,
2014
:
Offline estimates and tuning of mesospheric gravity-wave forcing using Met Office analyses
.
Quart. J. Roy. Meteor. Soc.
,
140
,
1025
1038
, https://doi.org/10.1002/qj.2168.
Lott
,
F.
, and
L.
Guez
,
2013
:
A stochastic parameterization of the gravity waves due to convection and its impact on the equatorial stratosphere
.
J. Geophys. Res.
,
118
,
8897
8909
, https://doi.org/10.1002/JGRD.50705.
Manzini
,
E.
, and
N. A.
McFarlane
,
1998
:
The effect of varying the source spectrum of a gravity wave parameterization in a middle atmosphere general circulation model
.
J. Geophys. Res.
,
103
,
31 523
31 539
, https://doi.org/10.1029/98JD02274.
McCormack
,
J. P.
,
S. D.
Eckermann
, and
T. F.
Hogan
,
2015
:
Generation of a quasi-biennial oscillation in an NWP model using a stochastic gravity wave drag parameterization
.
Mon. Wea. Rev.
,
143
,
2121
2147
, https://doi.org/10.1175/MWR-D-14-00208.1.
McCormack
,
J. P.
, and Coauthors
,
2017
:
Comparison of mesospheric winds from a high-altitude meteorological analysis system and meteor radar observations during boreal winters of 2009–2010 and 2012–2013
.
J. Atmos. Sol.-Terr. Phys.
,
154
,
132
166
, https://doi.org/10.1016/j.jastp.2016.12.007.
McLandress
,
C.
,
1998
:
On the importance of gravity waves in the middle atmosphere and their parameterization in general circulation models
.
J. Atmos. Sol.-Terr. Phys.
,
60
,
1357
1383
, https://doi.org/10.1016/S1364-6826(98)00061-3.
Osprey
,
S. M.
,
N.
Butchart
,
J. R.
Knight
,
A. A.
Scaife
,
K.
Hamilton
,
J. A.
Anstey
,
V.
Schenzinger
, and
C.
Zhang
,
2016
:
An unexpected disruption of the atmospheric quasi-biennial oscillation
.
Science
,
353
,
1424
1427
, https://doi.org/10.1126/science.aah4156.
Peña-Ortiz
,
C.
,
H.
Schmidt
,
M. A.
Giorgetta
, and
M.
Keller
,
2010
:
QBO modulation of the semiannual oscillation in MAECHAM5 and HAMMONIA
.
J. Geophys. Res.
,
115
,
D21106
, https://doi.org/10.1029/2010JD013898.
Pulido
,
M.
,
S.
Polavarapu
,
T.
Shepherd
, and
J.
Thuburn
,
2012
:
Estimation of optimal gravity wave parameters for climate models using data assimilation
.
Quart. J. Roy. Meteor. Soc.
,
138
,
298
309
, https://doi.org/10.1002/qj.932.
Richter
,
J. H.
,
F.
Sassi
, and
R. R.
Garcia
,
2010
:
Toward a physically based gravity wave source parameterization in a general circulation model
.
J. Atmos. Sci.
,
67
,
136
156
, https://doi.org/10.1175/2009JAS3112.1.
Ridout
,
J. A.
,
Y.
Jin
, and
C.-S.
Lou
,
2005
:
A cloud-base quasi-balance constraint for parameterized convection: Application to the Kain–Fritsch cumulus scheme
.
Mon. Wea. Rev.
,
133
,
3315
3334
, https://doi.org/10.1175/MWR3034.1.
Scaife
,
A. A.
,
N.
Butchart
,
C. D.
Warner
,
D.
Stainforth
,
W.
Norton
, and
J.
Austin
,
2000
:
Realistic quasi-biennial oscillations in simulation of the global climate
.
Geophys. Res. Lett.
,
27
,
3481
3484
, https://doi.org/10.1029/2000GL011625.
Scaife
,
A. A.
, and Coauthors
,
2014
:
Predictability of the quasi-biennial oscillation and its northern winter teleconnection on seasonal to decadal timescales
.
Geophys. Res. Lett.
,
41
,
1752
1758
, https://doi.org/10.1002/2013GL059160.
Schirber
,
S.
,
E.
Manzini
, and
M. J.
Alexander
,
2014
:
A convection-based gravity wave parameterization in a general-circulation model: Implementation and improvements on the QBO
.
J. Adv. Model. Earth Syst.
,
6
,
264
279
, https://doi.org/10.1002/2013MS000286.
Scinocca
,
J. F.
,
N. A.
McFarlane
,
M.
Lazare
,
J.
Li
, and
D.
Plummer
,
2008
:
The CCCma third generation AGCM and its extension into the middle atmosphere
.
Atmos. Chem. Phys.
,
8
,
7055
7074
, https://doi.org/10.5194/acp-8-7055-2008.
Serva
,
F.
,
C.
Cagnazzo
,
A.
Riccio
, and
E.
Manzini
,
2018
:
Impact of a stochastic nonorographic gravity wave parameterization on the stratospheric dynamics of a general circulation model
.
J. Adv. Model. Earth Syst.
,
10
,
2147
2162
, https://doi.org/10.1029/2018MS001297.
Shutts
,
G.
, and
S.
Vosper
,
2011
:
Stratospheric gravity waves revealed in NWP forecast models
.
Quart. J. Roy. Meteor. Soc.
,
137
,
303
317
, https://doi.org/10.1002/qj.763.
Skamarock
,
W. C.
,
2004
:
Evaluating mesoscale NWP models using kinetic energy spectra
.
Mon. Wea. Rev.
,
132
,
3019
3032
, https://doi.org/10.1175/MWR2830.1.
Tandeo
,
P.
,
M.
Pulido
, and
F.
Lott
,
2015
:
Offline estimation of subgrid-scale orographic parameters using EnKF and maximum likelihood error covariance estimates
.
Quart. J. Roy. Meteor. Soc.
,
141
,
383
395
, https://doi.org/10.1002/qj.2357.
Tsuda
,
T.
,
M.
Nishia
,
C.
Rocken
, and
R. H.
Ware
,
2000
:
A global morphology of gravity wave activity in the stratosphere revealed by the GPS occultation data (GPS/MET)
.
J. Geophys. Res.
,
105
,
7257
7273
, https://doi.org/10.1029/1999JD901005.
Vosper
,
S. B.
,
A. R.
Brown
, and
S.
Webster
,
2016
:
Orographic drag on islands in the NWP mountain grey zone
.
Quart. J. Roy. Meteor. Soc.
,
142
,
3128
3137
, https://doi.org/10.1002/qj.2894.
Xue
,
X.-H.
,
H.-L.
Liu
, and
X.-K.
Dou
,
2012
:
Parameterization of the inertial gravity waves and generation of the quasi-biennial oscillation
.
J. Geophys. Res.
,
117
,
D06103
, https://doi.org/10.1029/2011JD016778.
Yang
,
G.-Y.
,
B. J.
Hoskins
, and
J. M.
Slingo
,
2011
:
Equatorial waves in opposite QBO phases
.
J. Atmos. Sci.
,
68
,
839
862
, https://doi.org/10.1175/2010JAS3514.1.
Yao
,
W.
, and
C.
Jablonowski
,
2015
:
Idealized quasi-biennial oscillations in an ensemble of dry GCM dynamical cores
.
J. Atmos. Sci.
,
72
,
2201
2226
, https://doi.org/10.1175/JAS-D-14-0236.1.
Yu
,
C.
,
X.
Xue
,
J.
Wu
,
T.
Chen
, and
H.
Li
,
2017
:
Sensitivity of the quasi-biennial oscillation simulated in WACCM to the phase speed spectrum and the settings in an inertial gravity wave parameterization
.
J. Adv. Model. Earth Syst.
,
9
,
389
403
, https://doi.org/10.1002/2016MS000824.
Zhang
,
Y.
,
J.
Xiong
,
L.
Liu
, and
W.
Wan
,
2012
:
A global morphology of gravity wave activity in the stratosphere revealed by the 8-year SABER/TIMED data
.
J. Geophys. Res.
,
117
,
D21101
, https://doi.org/10.1029/2012JD017676.

Footnotes

For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).