## 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).

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 *c*_{j} in the range *u*_{src} ± *c*_{r}, where *u*_{src} is the zonal wind at the launch level (model pressure *p*_{src}) and *c*_{r} is called the phase speed range. By design, the scheme does not parameterize waves with source-level intrinsic phase speeds |*c*_{j} − *u*_{src}| greater than |*c*_{r}|. The source-level momentum flux of wave *j* is a Gaussian function of intrinsic phase speed:

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

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

where *k* is zonal wavenumber, *ρ* is density, *u* is the background zonal wind at that level, Fr_{c} 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:

where *g* is gravitational acceleration and *ϵ* is a parameterized efficiency (or intermittency) of wave breaking. The acceleration *a*_{j} is given the same sign as *u* − *c*_{j}. Note that a parameterized wave is totally absorbed at its critical level where *c*_{j} = *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 *a*_{j} 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 *u*_{src} ± *c*_{r}. 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 *p*_{src} (default 500 hPa), phase speed Gaussian width *c*_{0} (default 30 m s^{−1}), phase speed range *c*_{r} (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 *a*_{j} 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 (|*c*_{j} − *u*_{src}| ~ 30 m s^{−1}) waves will see a significant absolute change in source momentum flux, while waves with phase speeds nearer *u*_{src} 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 *a*_{j} 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}, *p*_{src}, *c*_{0}, *c*_{r}, 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\xaf$. 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\xaf$ profile over the hindcast period Δ*t* of $\Delta u\xaf$, we can then express the model response as a simple linear matrix $M$ governing the perturbation equation:

where **x**′ is the vector of perturbations of the GWD source parameters from their default values and $\Delta u\xafdef$ is the change in $u\xaf$ 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 *t*_{0}, $\Delta u\xaf\u2032$ simplifies to $(u\xaf\u2212u\xafdef)|t=t0+\Delta t=u\xaf\u2032$, 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 $\Delta u\xaf$ 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 *A* − *F* error) we first replace $u\xaf|t=t0+\Delta t$ above with the state we require, $u\xafA$, 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\xafA$ 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 $\Vert Mx\u2032\u2212\u2061(u\xafA\u2212u\xafdef)|t=t0+\Delta t\Vert $.

One implicit assumption of such an approach is that all of the *A* − *F* 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 $\Delta u\xaf$ 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:

where each column in $X$ is the vector of perturbations of the GWD source parameters in the selected subset for ensemble member *e* ($xe\u2032$), and the corresponding column in $U$ is the deviation of the hindcast profile of ensemble member *e* ($u\xafe$) from the default hindcast ($u\xafdef$) at time *t* = *t*_{0} + Δ*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\xafe\u2032$, with the perturbation profile predicted using the linear model $M$ and the ensemble member parameters $xe\u2032$ [i.e., $u\xafpred\u2032=Mxe\u2032$ as in Eq. (4)]. The squared correlation coefficient (*R*^{2}) between $u\xafe\u2032$ and $u\xafpred\u2032$ will be used to assess the accuracy of the linearization approximation, with high *R*^{2} values denoting strong linearity. Low *R*^{2} values may indicate that interaction terms, which we assumed did not exist when we wrote Eq. (4), contribute significantly to $u\xaf\u2032$. Poor prediction skill in the region of interest may also promote low *R*^{2} 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 $\Delta u\xaf$ 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 $\Delta u\xafdef$ 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 3–7 detail results from the ensemble of equatorial $u\xaf$ 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\xaf$ 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\xafe\u2032$) with the predicted response $u\xafpred\u2032=Mxe\u2032$ as in Eq. (4) over all analysis levels between 70 and 1 hPa (16 levels).

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\xaf\u2032$. 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 *c*_{0} and *τ*_{src} of about 24 m s^{−1} and 3.0 mPa, respectively. The *R*^{2} values are calculated to be 0.941 and 0.987, which again verifies that a single-predictor model using *c*_{0} or *τ*_{src} is able to accurately reproduce hindcast profiles of $u\xaf$ using the perturbed parameter value and our tangent-linear matrix $M$. Note also that changes in *c*_{0} 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 *τ*_{src}*c*_{0}. Equivalent changes to *τ*_{src} and *c*_{0} 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 *c*_{0} scales relative to the parameterized wave’s intrinsic phase speed. Thus, we expect that changes in *c*_{0} (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 *c*_{r} behaves uniquely among the five source parameters (Fig. 5), in that it appears to exhibit linearity only within a small window between *c*_{r} = 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 *c*_{r} values. Consider the case of *c*_{r} < 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 *c*_{r} 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 *c*_{r} 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 *p*_{src} (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 *p*_{src} and *k* on the hindcast $u\xaf$ profiles are too minimal for a realistic value of either parameter to provide a measurable reduction in *A* − *F* error. Note that this conclusion refers only to errors in $u\xaf$ 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}, *p*_{src}, *c*_{0}, *c*_{r}, 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. 3–7) 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\xaf\u2032$, 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., $\Delta u\xaf$) 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., $\Delta u\xaf=u\xafoptimal\u2212u\xafdef=Mx*$ as in Eq. (4)] is shown in red, and the validating analysis $\Delta u\xaf$ 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 *p*_{src} 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\xafpred\u2032=Mxe\u2032$ and the actual hindcast ensemble member perturbation profile $u\xafe\u2032$ has an *R*^{2} value of 0.554, indicating only moderate confidence that the model could closely reproduce the hindcast wind profile $u\xafe$ of ensemble member *e* when estimated using only $M$ and $xe\u2032$ as in Eq. (4). Since no pair of predictors is uncorrelated and nonlinear terms are likely to affect $u\xaf\u2032$ significantly, a model using only a single parameter (i.e., Figs. 3–7) is the best choice for the optimization procedure when applied to equatorial stratospheric $u\xaf$.

### d. Methodology validation

Large *R*^{2} values calculated using our linear matrix model [Eq. (4)] with *c*_{0} or *τ*_{src} as the lone predictor (Figs. 3 and 4) show that we are able to reproduce ensemble member hindcast $u\xaf$ using $M$ and $xe\u2032$. 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\xaf$ profile estimated by evaluating Eq. (4) using *τ**** and the $M$ for that initial date, while the green curve shows the corresponding $u\xaf$ 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.

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\xaf\u2032$ will use phase speed width *c*_{0} or source strength *τ*_{src} as a single predictor. We have validated the optimization methodology by demonstrating that our assumption that $u\xaf\u2032$ varies linearly with **x**′ is justifiable, that the linear model $M$ can faithfully reproduce $u\xaf'$ 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 {*c*_{0}, *τ*_{src}}, we observe collinearity of the response profiles and an *R*^{2} 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 *c*_{0} 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\xaf$ 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 *c*_{0} 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 *A* − *F* 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).

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 *A* − *F* 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 *R*^{2} 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 *c*_{0} 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 *c*_{0} 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 *R*^{2} decreased as the lead time increased (see Fig. 11c). However, the lower stratosphere fit (green curve) continues to exhibit a relatively high *R*^{2} 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.

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 *c*_{0} 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, *c*_{0} is more effective at increasing forecast skill than *τ*_{src} at both subseasonal and seasonal time scales; however, it also suffers from lower *R*^{2}.

### 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 *R*^{2} 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 *A* − *F* 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 (*R*^{2} > 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 *R*^{2} > 0.6 for any month at any longer lead time.

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\xaf$ 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.

## 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).

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 *A* − *F* 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 *c*_{0} 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 *c*_{0} 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 *c*_{0} 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

## Footnotes

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