Abstract

The climatology and future changes of the Brewer–Dobson circulation (BDC) in three climate change scenarios are studied using the latest version of the Whole Atmosphere Community Climate Model (WACCM4), which is fully coupled to an ocean model. The results show an acceleration in both the shallow and deep branches of circulation in response to increasing greenhouse gases (GHGs) together with an upward displacement of the tropical upwelling in the deep branch near the stratopause. The downward control principle reveals that different waves are involved in forcing the acceleration of the upper and lower branches. Climatological-mean tropical upwelling in both the lower and upper stratosphere is dominated by explicitly resolved, planetary-scale waves. Trends in the tropical upwelling in the lower stratosphere are mainly attributed to explicitly resolved, planetary-scale waves. However, in the upper stratosphere, despite the fact that resolved waves control the forcing of the climatological upwelling, their contribution to the long-term trend diminishes with increasing GHGs, while the role of gravity waves associated with fronts increases and becomes dominant in the model scenario with the largest GHG increases. The intensification and upward displacement of the subtropical tropospheric jets due to climate change leads to filtering of the westerly part of the frontal gravity wave spectrum, leaving the easterly components to reach the upper stratosphere and force the changes in the circulation there.

1. Introduction

The stratospheric meridional overturning circulation, also known as the Brewer–Dobson circulation (BDC), is responsible for mean mass transport of chemical tracers from the equator to the poles (Birner and Bönisch 2011). Several modeling studies have shown an acceleration of the BDC with increasing greenhouse gas (GHG) concentrations (e.g., Butchart et al. 2006), which led to growing interest in this circulation. Changes in the BDC can modify the distribution of ozone and other chemical compounds, resulting not only in a dynamical alteration of the stratosphere and stratosphere–troposphere interactions but also to changes in the tropospheric radiative balance.

Recent studies have clearly differentiated two branches of the stratospheric mean meridional circulation: a shallow branch, located in the lowermost stratosphere, with upwelling in the tropics and downwelling in the subtropics and middle latitudes, and a deep branch with maximum upwelling in the tropical upper stratosphere and downwelling in the middle and high latitudes throughout the entire height of the stratosphere. Birner and Bönisch (2011) distinguish between shallow and deep branches on the basis of the stratospheric transit times in each branch, with the result that the shallow branch extends to about 50 hPa and the deep branch is located above that altitude. Here we follow Birner and Bönisch and adopt 50 hPa as the approximate level separating shallow and deep branches.

Recently, Lin and Fu (2013) looked at the climatology and changes in all branches of the BDC using chemistry–climate models from the second Chemistry–Climate Model Validation (CCMVal2) activity. These authors define a transition layer between 100 and 70 hPa, with the shallow branch extending up to 30 hPa, and the deep branch beginning above 30 hPa. They found a larger acceleration in the shallow branch than in the deep branch. Similarly, Hardiman et al. (2014) investigated future changes in the Brewer–Dobson circulation in the entire stratosphere in high-top phase 5 of the Coupled Model Intercomparison Project (CMIP5) models (with model lids above 1 hPa). They found differences in the width of the changes in tropical upwelling above and below 20 hPa. These differences suggest that different dynamical behaviors may be related to different driving mechanisms.

The stratospheric mean meridional circulation is driven by wave dissipation (Holton et al. 1995). The downward control principle (DCP; Haynes et al. 1991) makes it possible to identify quantitatively the waves that are responsible for driving this circulation. This principle establishes that, in steady state, the mean vertical velocity at any vertical level is controlled only by the density-weighted integral of the wave forcing acting above that level. Different studies using both reanalysis and general circulation models suggest that the BDC in the lower stratosphere is mainly forced by large-scale waves [e.g., Fig. 4.10 of Eyring et al. (2010)], although there is considerable intermodel variability. In particular, in CCMVal2 models, more than 70% of the climatological tropical upwelling at 70 hPa was due to waves that are explicitly resolved by the models, 21.1% to parameterized (unresolved) orographic gravity wave drag, and 7.1% to parameterized nonorographic gravity wave drag (Butchart et al. 2011). These percentages agree well with results from the 1989–2009 Interim European Centre for Medium-Range Weather Forecasts (ECMWF) Re-Analysis (ERA-Interim), although in the latter only 4% of the tropical upwelling at 70 hPa was attributed to parameterized wave drag (nonorographic gravity waves are not included in ERA-Interim; Seviour et al. 2012).

Regarding changes in tropical upwelling associated with increasing greenhouse gases, there is a general consensus that changes in resolved wave drag are the main driver of the acceleration in the lower stratosphere [e.g., Fig. 4.11 of Eyring et al. (2010)]. The multimodel analysis of Butchart et al. (2010) revealed contributions at 70 hPa of up to 67% from these waves and 30% from orographic gravity waves. These results agree with other individual model studies, such as those by Garcia and Randel (2008; hereafter GR08), McLandress and Shepherd (2009), Garny et al. (2011), and Shepherd and McLandress (2011). GR08 pointed out that, as the subtropical tropospheric jets intensify in response to tropospheric warming driven by GHG, the contribution of gravity waves increases above 70 hPa in the middle stratosphere. This altitude dependence of the different wave contributions to the upwelling highlights the importance of analyzing the two branches of the circulation separately.

Several studies have attempted to determine whether changes in wave forcing that drive changes in the circulation are due to differences in wave transmission or in wave excitation. In particular, Calvo and Garcia (2009) found that changes in the transmission of explicitly resolved waves were the main mechanism in simulations of the last half of the twentieth century. This was in contrast with simulations of the twenty-first century, where changes in the excitation of resolved waves, related to anomalous convection, dominated over changes in wave transmission. Olsen et al. (2007) showed changes in tropospheric wave propagation due to warmer sea surface temperatures (SSTs) in the tropics. Furthermore, Deckert and Dameris (2008) demonstrated that higher tropical SSTs enhance tropical deep convection, amplifying the generation of tropical waves, as well as poleward-propagating waves that can intensify the BDC in the extratropics. Finally, Garny et al. (2011) showed that both wave propagation and wave generation are affected by increasing SSTs as a result of climate change and, thus, both mechanisms can generate changes in the BDC. In particular, Shepherd and McLandress (2011) found that changes in wave drag in their model were largely explained in terms of changes in the location of the critical layers within the subtropical lower stratosphere.

The fact that SSTs represent an important atmosphere–ocean feedback and that its explicit simulation in climate models is complex has been noted in many previous studies (e.g., Braesicke and Pyle 2004; Kosaka and Xie 2013). Up to now, most studies that investigated changes in the BDC with general circulation models that include a well-resolved middle atmosphere used prescribed SSTs from observations for past simulations or SSTs from low-top atmosphere–ocean models to understand future changes associated with GHGs. This lack of interactive atmosphere–ocean feedbacks adds uncertainty to the understanding of the mechanisms involved in the BDC changes. In the present study, we use a fully coupled atmosphere–ocean climate model with its top in the thermosphere to analyze the BDC and its forcing in the entire stratosphere. Special attention is paid to the less-investigated deep branch of the circulation and to future changes under different climate change scenarios.

This paper is organized as follows: Section 2 presents the model and simulations used here together with the methodology. Results concerning the BDC climatology are shown in section 3, while section 4 presents future trends in the circulation. Section 5 summarizes the main results.

2. Simulations and method

The latest version of the Whole Atmosphere Community Climate Model (WACCM4) is part of the Community Earth System Model, version 1 (CESM1), developed by the U.S. National Center for Atmospheric Research (NCAR). This improved version of WACCM integrates a four-component coupled system, including atmosphere, ocean, land, and sea ice models. The ocean and sea ice components are described in Holland et al. (2012) and Danabasoglu et al. (2012). The atmospheric component, which includes fully coupled chemistry, has a longitude–latitude resolution of 2.5° × 1.9° and 66 vertical levels up to about 140 km. The vertical resolution is finer than 1.5 km in the first 25 km, with a gradual increase up to 2 km at the stratopause, 3.5 km in the mesosphere, and half of the local scale height beyond the mesopause. The model includes heating from volcanic aerosols (Tilmes et al. 2009), solar variability is specified from the model of Lean et al. (2005), and a quasi-biennial oscillation is imposed by relaxing the winds to observations in the tropics (Matthes et al. 2010). The parameterization of orographic gravity waves is related to surface roughness, and nonorographic gravity waves are parameterized based on the diagnosis of frontogenesis and the occurrence of deep convection in the model. These parameterizations are described in detail by Richter et al. (2010). Additional information about the model can be found in Marsh et al. (2013).

Three historical runs were performed from 1960 to 2005, each with slightly different initial conditions, and observed concentrations of GHGs at the surface. Three future simulations from 2005 to 2100 were run based on the representative concentration pathways (RCPs) adopted by the Intergovernmental Panel on Climate Change (IPCC) for its Fifth Assessment Report. RCP2.6 is the most optimistic scenario and prescribes a decrease in CO2 starting in 2050, reaching 400 ppmv in 2100. In RCP4.5, the CO2 concentration increases from 400 ppmv to slightly less than 550 ppmv in 2100. The RCP8.5 simulation represents a CO2 doubling from 2005 to 2080, reaching values of 950 ppmv in 2100. Mixing ratios of ozone depleting substances (ODS) in terms of equivalent effective stratospheric chlorine (EESC) decrease from maximum values of approximately 4 ppbv in 2005 to approximately 1.7 ppbv in 2100.

To characterize the mean meridional circulation of the stratosphere, we examine the vertical component of the transformed Eulerian-mean (TEM) residual circulation (Andrews et al. 1987). We will refer to upwelling when is positive and to downwelling when is negative. Possible changes in the future are assessed by calculating linear trends derived from deseasonalized monthly-mean time series. The significance of trends at the 95% confidence level is evaluated using a Monte Carlo test. Differences between climatological values from the first and last 20 years of the simulation are also shown, as needed. Significant differences are assessed at the 95% confidence level by a Student’s t test.

3. Climatology of the BDC in WACCM4

a. Circulation

Figure 1a shows the latitude–altitude cross section of the annual-mean for the ensemble mean of the three historical simulations (1960–2005). This figure clearly distinguishes the two different circulation branches discussed in the introduction. The shallow branch, in the lower stratosphere, is principally responsible for troposphere–stratosphere mass exchange and consists of two cells with rising motion in the tropics and sinking in the subtropics and middle latitudes (Holton 1990). Tropical upwelling in the shallow branch is stronger in December–February (DJF) than in June–August (JJA). In the deep branch, the largest tropical upwelling occurs in the upper stratosphere around 45–50 km. The associated downwelling extends from midlatitudes to the poles, reaching a maximum at polar latitudes in the upper stratosphere and extending down toward the polar troposphere. Similar to the shallow branch, the deep branch is equatorially centered when the annual average is computed. However, the seasonal averages show an intensification of the winter cell and a shift of the tropical upwelling toward the summer hemisphere, also stronger in DJF (Seviour et al. 2012). The deep branch cell in the summer hemisphere is very weak below 40 km, in contrast with the behavior of the shallow branch (Figs. 1b,c). Note that the summer to winter hemisphere branch of the TEM circulation starts above 65–70 km in WACCM4, beyond the upper boundary of the figures shown here.

Fig. 1.

Climatology of the TEM (, ) vector circulation (arrows) and magnitude of (contours) as functions of latitude and log-pressure altitude zp for the ensemble mean of three simulations of the historical period 1960–2005. (a) Annual-mean climatology, (b) DJF mean, and (c) JJA mean. For , variable contours are shown for 0, ±0.1, ±0.2 mm s−1, every 0.2 mm s−1 up to ±1 mm s−1, and every 1 mm s−1 thereafter. Orange shading and solid lines represent positive values; blue shading and dashed lines represent negative values. The arrows on the lower-left-hand corner of each panel indicate the scale of the vector circulation.

Fig. 1.

Climatology of the TEM (, ) vector circulation (arrows) and magnitude of (contours) as functions of latitude and log-pressure altitude zp for the ensemble mean of three simulations of the historical period 1960–2005. (a) Annual-mean climatology, (b) DJF mean, and (c) JJA mean. For , variable contours are shown for 0, ±0.1, ±0.2 mm s−1, every 0.2 mm s−1 up to ±1 mm s−1, and every 1 mm s−1 thereafter. Orange shading and solid lines represent positive values; blue shading and dashed lines represent negative values. The arrows on the lower-left-hand corner of each panel indicate the scale of the vector circulation.

A different look into the climatological seasonal behavior of the BDC is presented in Fig. 2. It shows the latitudinal distribution of the ensemble-mean, annual-mean as a function of time and latitude at 1 hPa (Fig. 2a) and 100 hPa (Fig. 2b). These are the levels where the tropical upwelling in WACCM4 is largest in the deep and shallow branches, respectively (Fig. 1a). The seasonal behavior of is similar in both the lower and the upper stratosphere, although the seasonal cycle is much more pronounced in the upper stratosphere (Fig. 2a). Thus, at 1 hPa the upwelling can comprise the entire winter hemisphere while downwelling occurs in the summer hemisphere from about 15° latitude to the pole. The magnitude of the tropical upwelling is very similar throughout the year, slightly larger in boreal winter months. In the Southern Hemisphere (SH), the largest downwelling occurs at middle latitudes and not in the polar region. In the lower stratosphere, upwelling is confined to lower latitudes, maximizing in the summer hemisphere in the subtropics. Downwelling is simulated from the subtropics to the polar latitudes, and it is larger in the winter half of the year. Note, however, that downwelling at polar latitudes in the lower stratosphere comes mainly from the deep branch (e.g., Birner and Bönisch 2011). The seasonal distribution of at 100 hPa can be compared with results from the ERA-Interim, as shown by Seviour et al. (2012, their Fig. 4c). Although the period is not exactly the same (1989–2009 for the reanalysis versus 1960–2005 for the model) the agreement is excellent, not only regarding the spatial distribution of the regions of upwelling and downwelling, but also their magnitude. This good agreement lends confidence to our results and supports the use of WACCM4 to investigate future changes in the BDC.

Fig. 2.

Ensemble-mean seasonal climatology of (1960–2005) at (a) 1 and (b) 100 hPa. Contour intervals are 0.2 mm s−1 up to 1 mm s−1 and every 1 mm s−1 thereafter. Orange shading and solid lines represent positive values; blue shading and dashed lines represent negative values.

Fig. 2.

Ensemble-mean seasonal climatology of (1960–2005) at (a) 1 and (b) 100 hPa. Contour intervals are 0.2 mm s−1 up to 1 mm s−1 and every 1 mm s−1 thereafter. Orange shading and solid lines represent positive values; blue shading and dashed lines represent negative values.

b. Wave forcing

In WACCM4, the wave drag from resolved and parameterized waves is output separately. Within the gravity wave parameterization, orographic gravity waves and nonorographic gravity waves of frontal and convective origin are also distinguished. Thus, it is possible to use the DCP to quantify the different contribution of these waves to the tropical upwelling. The DCP estimate of averaged between latitudes and is given by

 
formula

where ρ is the density, z is log-pressure altitude, is the zonal-mean angular momentum, is Earth’s rate of rotation, and a is the radius of Earth. The terms in the integrand include the Eliassen–Palm flux divergence due to resolved waves, · F, and the parameterized wave drag X due to orographic gravity waves and gravity waves associated with convective and frontal sources. We average (1) over the turn-around latitudes (i.e., and are taken to be the latitudes where changes sign on each pressure level; Rosenlof 1995). In what follows, the individual parameterized wave drag components will be denoted by OGWD, CGWD, and FGWD, which stand for orographic, convective, and frontal gravity wave drag, respectively. This methodology is similar to that used in GR08 except that, for each altitude, we evaluate the integrand in (1) between the turn-around latitudes instead of using fixed latitudes as in GR08. Note that, in order for the DCP to be applicable, the turn-around latitudes where the integrand is evaluated need to be approximately parallel to lines of constant angular momentum, and thus (1) is not valid in the deep tropics where lines of constant m form closed contours.

Figure 3 shows computed from (1) for the ensemble mean of the three historical simulations. The latitude limits used in the calculation correspond to the annual-mean climatological turn-around latitudes for 1960–2005. These latitude limits are far enough from the equator that application of (1) is justified. The estimated value for computed from the DCP (dashed line) and its actual value output from the model (black line) are in very good agreement everywhere except near 48 km, although the difference there is only about 15%. The tropical upwelling clearly shows the two maxima, in the lowermost and upper stratosphere, corresponding to the shallow and deep branches of the circulation. The largest values (about 1 mm s−1) occur at about 48 km. Throughout the entire stratosphere, resolved waves are the main driver of the climatological tropical upwelling. In the lower stratosphere, at 70 hPa (~18.6 km), their contribution is about 70% of the total, versus 30% for orographic gravity waves. The contributions of convective gravity waves (CGWs) or frontal gravity waves (FGWs) are negligible in the lower stratosphere. These results are consistent with other modeling studies and reanalysis data (e.g., Butchart et al. 2011; Seviour et al. 2012).

Fig. 3.

Tropical average of calculated from the DCP, showing the contributions due to different waves in the 1960–2005 ensemble-mean simulation. The black curve includes the total wave drag, the red curve is the resolved wave drag, the blue curve is FGWD, the green curve is OGWD, and the yellow curve is CGWD. The dashed curve is the TEM vertical velocity computed directly from the model output. Similar results hold for all other simulations discussed in this paper.

Fig. 3.

Tropical average of calculated from the DCP, showing the contributions due to different waves in the 1960–2005 ensemble-mean simulation. The black curve includes the total wave drag, the red curve is the resolved wave drag, the blue curve is FGWD, the green curve is OGWD, and the yellow curve is CGWD. The dashed curve is the TEM vertical velocity computed directly from the model output. Similar results hold for all other simulations discussed in this paper.

In the upper stratosphere, at about 48 km (where the largest upwelling occurs), most of the upwelling is attributable to resolved waves. In fact, upwelling due to resolved waves is actually larger than the total upwelling near the stratopause because FGWs produce downwelling at this level. Both OGWD and CGWD make a small contribution to upwelling near the stratopause—about 10% of total upwelling. Note that in WACCM4 the contribution of CGWs to the mean tropical upwelling is close to zero throughout the stratosphere. The wave contribution to the climatological-mean tropical upwelling in future scenarios (not shown) is similar to the one shown in Fig. 3 for the recent past.

4. Future trends in the BDC due to increases in GHG

a. Circulation trends

Long-term trends, 2005–2100, in for three future scenarios are shown in Fig. 4, with the trend in the vector TEM circulation (, ) superimposed (arrows). The spatial configuration of the trends can be compared with that of the 1960–2005 climatology shown in Fig. 1a, as the climatologies for the future simulations showed very similar spatial patterns. Figure 4 shows an intensification of both branches of the circulation with increasing GHG concentrations. The intensification is very weak in RCP2.6, where only the change in the upper branch near the stratopause is significant. In RCP4.5, there is a moderate, statistically significant acceleration of the circulation in the regions of maximum upwelling and downwelling in the upper stratosphere and lowermost stratosphere. RCP8.5 shows the largest trends in the TEM circulation, both in the deep and the shallow branches. The upper-stratosphere trend implies an increase in over the period of analysis of approximately 50% of the climatological-mean value. Note that the regions of positive upwelling trends in the lower stratosphere and in the vicinity of the stratopause are narrower than the climatological upwelling (1960–2005) at those locations (red lines), indicating that the upwelling becomes more concentrated around the equator. This is in agreement with results (for the lower stratosphere) obtained from an ensemble of several high-top CMIP5 models (Hardiman et al. 2014). A common feature in all our simulations in the SH is the positive trend in at high latitudes, implying a deceleration of the climatological downwelling over the southern polar cap (cf. Fig. 1). Comparison of the trends in the RCP4.5 scenario (Fig. 4b) with trends derived from a pair simulations (not shown), where either the mixing ratios of ODS or those of GHG were held constant at 2001 values, revealed that the deceleration of downwelling over the southern polar cap is related to ozone recovery. This has also been pointed out in other, recent studies (e.g., McLandress and Shepherd 2009).

Fig. 4.

Trends in the TEM (, ) vector circulation (arrows) and (contours) as functions of latitude and zp over the period 2005–2100 for (a) RCP2.6, (b) RCP4.5, and (c) RCP8.5. Contours are 0, ±0.01, ±0.02 mm s−1 decade−1, and every ±0.02 mm s−1 decade−1 thereafter. Stippling indicates where trends in are not significant at the 2σ level. Orange shading and solid lines represent positive values; blue shading and dashed lines represent negative values. Red lines indicate the climatological (1960–2005) turn-around latitudes. The arrow on the lower-left-hand corner of each panel indicates the scale of the vector trend.

Fig. 4.

Trends in the TEM (, ) vector circulation (arrows) and (contours) as functions of latitude and zp over the period 2005–2100 for (a) RCP2.6, (b) RCP4.5, and (c) RCP8.5. Contours are 0, ±0.01, ±0.02 mm s−1 decade−1, and every ±0.02 mm s−1 decade−1 thereafter. Stippling indicates where trends in are not significant at the 2σ level. Orange shading and solid lines represent positive values; blue shading and dashed lines represent negative values. Red lines indicate the climatological (1960–2005) turn-around latitudes. The arrow on the lower-left-hand corner of each panel indicates the scale of the vector trend.

In addition to the intensification of the BDC shown in Fig. 4, the largest trends in tropical upwelling in the upper stratosphere indicate an upward displacement of the region of upwelling with increasing GHG concentrations. Thus, for simulation RCP2.6, the maximum trend in upwelling is centered near 48 km, while for RCP8.5, it occurs near 52 km. To analyze this behavior in more detail, differences between the averages of the last and the first 20 years of the simulations were computed. The results are shown in Fig. 5 for RCP8.5, which is the scenario with the largest changes. Significant differences at the 95% level were computed with a Student’s t test. The difference pattern (Fig. 5b) in the upper stratosphere is centered at a higher altitude than the climatological upwelling. This upward displacement is larger in RCP8.5 than in RCP4.5, and it is not present in the RCP2.6 simulation (not shown), which indicates that certain changes only occur for the larger increases in GHGs concentrations. The behavior illustrated in Fig. 5 may be understood as an upward expansion of the deep branch tropical upwelling.

Fig. 5.

(a) averaged over the first 20 years of the RCP8.5 simulation (2005–25); contour interval is 0.5 mm s−1. (b) Difference in averaged over 2080–99 minus 2005–25. Variable contours are 0, ±0.05, ±0.1, ±0.2 mm s−1, and every 0.2 mm s−1 thereafter. Stippling indicates where trends in are not significant at the 2σ level according to a Student’s t test. Orange shading and solid lines represent positive values; blue shading and dashed lines represent negative values.

Fig. 5.

(a) averaged over the first 20 years of the RCP8.5 simulation (2005–25); contour interval is 0.5 mm s−1. (b) Difference in averaged over 2080–99 minus 2005–25. Variable contours are 0, ±0.05, ±0.1, ±0.2 mm s−1, and every 0.2 mm s−1 thereafter. Stippling indicates where trends in are not significant at the 2σ level according to a Student’s t test. Orange shading and solid lines represent positive values; blue shading and dashed lines represent negative values.

b. Trends in wave forcing

To investigate the contribution of the different waves to changes in the BDC, the latitude–altitude distribution of forcing is first examined. Figure 6 depicts the latitude–altitude cross sections of the climatologies and trends in the four types of wave forcing studied here [Eliassen–Palm (EP) flux divergence of resolved waves, OGWD, FGWD, and CGWD]. Results are shown for the RCP8.5 climate change scenario since it produces the largest changes. In the lowermost stratosphere, significant negative trends in EP flux divergence appear in the subtropics (Fig. 6b). These represent an intensification and upward extension of the climatological wave dissipation in this region (Fig. 6a). Significant trends in orographic wave forcing also appear in the lower stratosphere between 30° and 40°N, intensifying the climatological drag between 20 and 25 km and weakening it below 20 km (Figs. 6c,d). Nonorographic gravity waves do not dissipate in the lower stratosphere, in agreement with results from the DCP, which showed no contribution from FGWs and CGWs to the climatological tropical upwelling in this region. In addition, no changes in FWGD or CGWD were found in the future simulation in this region.

Fig. 6.

(left) RCP8.5 climatological annual means and (right) trends per decade for (a),(b) EP flux divergence, (c),(d) OGWD, (e),(f) FGWD, and (g),(h) CGWD. In the left panels, contours are every ±0.1 m s−1 day−1 up to ±0.5 m s−1 day−1 and then every ±1 m s−1 day−1. In the right panels, contours are every ±0.01 m s−1 day−1 decade−1 up to ±0.05 m s−1 day−1 decade−1 and then every ±0.1 m s−1 day−1 decade−1. Orange shading and solid lines represent positive values; blue shading and dashed lines represent negative values. Stippling indicates where trends are not significant at the 2σ level.

Fig. 6.

(left) RCP8.5 climatological annual means and (right) trends per decade for (a),(b) EP flux divergence, (c),(d) OGWD, (e),(f) FGWD, and (g),(h) CGWD. In the left panels, contours are every ±0.1 m s−1 day−1 up to ±0.5 m s−1 day−1 and then every ±1 m s−1 day−1. In the right panels, contours are every ±0.01 m s−1 day−1 decade−1 up to ±0.05 m s−1 day−1 decade−1 and then every ±0.1 m s−1 day−1 decade−1. Orange shading and solid lines represent positive values; blue shading and dashed lines represent negative values. Stippling indicates where trends are not significant at the 2σ level.

Significant trends in EP flux divergence are also found in the upper stratosphere. Significant negative trends enhance the climatological wave dissipation in the high latitudes of the NH and shift upward that one in the subtropics and midlatitudes of the SH, above 50 km. At lower altitudes, between 40 and 50 km, there are positive trends in EP flux divergence that weaken the climatological wave dissipation therein. Taken together, the trend patterns above and below 50 km may be interpreted as an upward shift in the region of the largest EP flux divergence with respect to climatology. There is also a significant positive trend in the SH at high latitudes, probably related to ozone recovery, as the abundance of ozone over Antarctica impacts resolved wave propagation (Li et al. 2008). Here, this is manifested as a positive EP flux divergence trend associated with a decrease in wave dissipation in the southern polar latitudes.

Significant negative trends in OGWD are simulated in the subtropical upper stratosphere of the NH, enhancing the climatological wave-induced acceleration. Significant negative (easterly) trends in FGWD are simulated in the subtropics and NH midlatitudes (Fig. 6f). They indicate a decrease of the climatological westerly FGWD in the subtropics and an enhancement of the climatological easterly component in the NH midlatitudes (Fig. 6e). These changes occur mainly during winter in each hemisphere, as explained below. The CGWD is the result of dissipation of gravity waves excited by deep convection in the tropics. These waves dissipate at much higher altitudes (mainly above 50 km; Fig. 6g) and, accordingly, they were not found to drive the climatological tropical upwelling in the stratosphere (Fig. 3). In terms of trends, weaker trends than in the case of FGWD are observed above 50 km in the tropical and subtropical regions, although they intensify during JJA (not shown).

To quantify the contribution of each of these waves to the trends in the climatological tropical upwelling, the DCP has been applied to the RCP4.5 and RCP8.5 simulations (Fig. 7) in a similar way as it was applied to the historical climatology. That is, we have applied (1) to each RCP scenario between the climatological annual-mean turn-around latitudes of the historical ensemble (1960–2005). No results are shown for RCP2.6 since no significant trends in were found in the tropical region in this case (Fig. 4a).

Fig. 7.

Trends in due to different waves between the turn-around latitudes of the climatological ensemble mean (1960–2005) for (a) RCP4.5 and (b) RCP8.5. Color lines are as in Fig. 3. Shadowing indicates the 2σ confidence limits.

Fig. 7.

Trends in due to different waves between the turn-around latitudes of the climatological ensemble mean (1960–2005) for (a) RCP4.5 and (b) RCP8.5. Color lines are as in Fig. 3. Shadowing indicates the 2σ confidence limits.

In the lower stratosphere, the net trends in tropical-mean upwelling in RCP4.5 and RCP8.5 seen in Fig. 7 are consistent with the results shown in Fig. 4. Trends in resolved wave drag and in OGWD are the only important driver of the enhanced upwelling of the shallow branch up to about 18 km (~70 hPa); above this level, orographic gravity waves (OGWs) play a nonnegligible role. In the lower stratosphere, the relative contribution of the different waves to the total trend is similar in both RCP scenarios, although the magnitude of the trends is larger in RCP8.5. Our results in this region are consistent with those obtained by GR08, who applied the DCP between fixed latitudes in a previous version of WACCM. In the upper stratosphere, RCP4.5 and RCP8.5 show significant trends in centered above 50 km (~1 hPa), above the location of maximum climatological (Fig. 1a), in agreement with the upward extension of the circulation pointed out above. Interestingly, the contribution of FGWs to the trend in increases dramatically from RCP4.5 to RCP8.5 compared to the increase driven by resolved EP flux divergence. Thus, while in RCP4.5, 36% of the total upwelling trend near 50 km is attributed to FGWs (and 64% to resolved waves), in RCP8.5 the FGW contribution increases up to 56% versus 30% for resolved waves, such that changes in FGW driving become the largest contributor to the trend in tropical upwelling in the upper stratosphere in the scenario with the largest increase in GHGs.

We have also applied the DCP in simulation RCP8.5 between the climatological turn-around latitudes and the pole in each hemisphere (Figs. 8a,b)—that is, over the latitude ranges where there is climatological downwelling. We note, however, that in these regions the trends in vertical velocity are not uniformly negative. This is particularly true of the SH, where poleward of 65°S the trend in vertical velocity is actually positive below about 45 km (cf. Fig. 4). In the NH (Fig. 8a), the mean downwelling trend between the northern turn-around latitude and the pole in the lower stratosphere is attributed to resolved waves and OGW; at 18 km, the percentages are 70% and 30%, respectively. In the upper stratosphere, 50% of the trend is attributed to FGWs while the other 50% is due to OGWs and resolved waves; above 50 km, the contribution of FGWs is dominant. CGWs contribute negligibly to the intensification of downwelling in the NH, as near zero trends were found everywhere.

Fig. 8.

Trends in attributable to different waves for simulation RCP8.5: (a) in the Northern Hemisphere, between the northern turn-around latitude and the North Pole; (b) in the Southern Hemisphere, between the southern turn-around latitude and the South Pole; and (c) between 65°S and the South Pole. Line colors are as in Fig. 3. Shadowing indicates the 2σ confidence limits.

Fig. 8.

Trends in attributable to different waves for simulation RCP8.5: (a) in the Northern Hemisphere, between the northern turn-around latitude and the North Pole; (b) in the Southern Hemisphere, between the southern turn-around latitude and the South Pole; and (c) between 65°S and the South Pole. Line colors are as in Fig. 3. Shadowing indicates the 2σ confidence limits.

In the SH, Fig. 8b shows that in the lower stratosphere, almost 100% of the downwelling trend is due to resolved waves. However, above 45 km, where the downwelling trend extends to polar latitudes, FGW, CGW, and resolved wave contributions increase reaching similar magnitudes at 53 km. Above this altitude, resolved waves become again the principal contributor. OGW were found to contribute negligibly to the SH downwelling trends. Finally, in order to determine the cause of the positive trend in vertical velocity over the southern polar cap at altitudes below 45 km (implying a weaker downwelling), we have applied the DCP between 65°S and the South Pole (Fig. 8c). The results show that resolved waves and OGW are responsible for this positive trend over the southern polar cap, with negligible contributions from other waves. This result is consistent with the pattern of trends shown for the EP flux divergence (Fig. 6b), which shows a decrease in resolved wave forcing below 50 km at high southern latitudes.

Changes in tropical upwelling are forced by changes in wave dissipation, which in turn are modulated by changes in the background winds. To understand the mechanisms that drive the changes in wave forcing with increasing GHGs, Fig. 9 shows the climatology (contours) and the trends (colors) of the zonal-mean zonal wind for the three RCP scenarios. There is a direct relationship between the intensification and upward shift of the subtropical jets and the increasing load of GHGs. The largest trends are present in the extreme scenario, RCP8.5 (Fig. 9c), consistent with the largest trends in meridional temperature gradient (not shown). In agreement with GR08 and Shepherd and McLandress (2011), the intensification of the subtropical jets enhances wave propagation and dissipation in the subtropical stratosphere (Fig. 7b). The negative trends in the SH present in all the RCP simulations reflect the impact of ozone recovery.

Fig. 9.

Trends (shaded) and climatology (black contours) of the zonal-mean zonal wind for three RCP scenarios: (a) RCP2.6, (b) RCP4.5, and (c) RCP8.5. Contour interval for the trend is ±0.2 m s−1 decade−1; for the climatology, the interval is ±10 m s−1. Orange shading and solid lines represent positive values; blue shading and dashed lines represent negative values. Stippling indicates where trends are not significant at the 2σ level. Note that the altitude range starts from the surface in these plots.

Fig. 9.

Trends (shaded) and climatology (black contours) of the zonal-mean zonal wind for three RCP scenarios: (a) RCP2.6, (b) RCP4.5, and (c) RCP8.5. Contour interval for the trend is ±0.2 m s−1 decade−1; for the climatology, the interval is ±10 m s−1. Orange shading and solid lines represent positive values; blue shading and dashed lines represent negative values. Stippling indicates where trends are not significant at the 2σ level. Note that the altitude range starts from the surface in these plots.

Changes in the zonal-mean zonal wind can affect the propagation and dissipation of gravity waves. We show here how these changes affect FGWs, which are the main contributor to the changes in tropical upwelling found in the upper stratosphere in the RCP8.5 simulation (see Fig. 6f). Figures 10a and 10b show vertical profiles of the climatological FGW eddy momentum flux associated with waves of easterly (red) and westerly (blue) phase velocity, together with the zonal-mean zonal wind (black) at 20°N and 20°S. These are approximately the latitude limits for the maximum trend in the deep branch upwelling, near the stratopause. The net FGW momentum flux (green curve) is zero at the source but quickly becomes positive below 20 km owing to the filtering of westerly components of the spectrum by the climatological subtropical jet (black curve). Above 20 km, where the climatological zonal-mean zonal wind turns negative, easterly components of the spectrum are attenuated sharply, and the net FGW momentum flux becomes negative. The negative net flux decreases above 50 km, giving rise to the climatological positive wave drag seen in Fig. 6e.

Fig. 10.

Annual-mean (a),(b) climatology and (c),(d) trends of the frontal gravity wave momentum flux (red for waves of easterly phase velocity, blue for waves of westerly phase velocity, and green for the net flux) at (left) 20°S and (right) 20°N. Also shown in (c),(d) are trends in zonal-mean zonal wind (m s−1 decade−1; black). Shadowing indicates the 2σ confidence limits.

Fig. 10.

Annual-mean (a),(b) climatology and (c),(d) trends of the frontal gravity wave momentum flux (red for waves of easterly phase velocity, blue for waves of westerly phase velocity, and green for the net flux) at (left) 20°S and (right) 20°N. Also shown in (c),(d) are trends in zonal-mean zonal wind (m s−1 decade−1; black). Shadowing indicates the 2σ confidence limits.

Trends in FGW momentum flux are shown in Figs. 10c and 10d. At the source of the waves, near 5 km, both westerly and easterly components of the spectrum become weaker with respect to climatology; that is, the tendencies have the opposite sign of the climatological values. This behavior is associated with a weakening of the zonal-mean zonal wind below 10 km. It is not altogether clear why this weakening of the tropospheric zonal-mean wind reduces the FGW momentum fluxes at their source, but a possibility is that the frequency of occurrence of fronts diminishes at this location; since the diagnosis of frontogenesis is the trigger for the FGWs in the model (Richter et al. 2010), infrequent frontogenesis will result in a weakening of the momentum flux associated with these waves. Some distance above the wave source, at 15–20 km, there is a sharp intensification of the momentum flux carried by the easterly components of the spectrum, which is most pronounced in the SH. This is related to the strengthening of the subtropical jets with increasing GHG concentrations (black curves). Above about 25 km, where the trends in the zonal-mean wind become negative, the flux carried by the easterly components is in turn attenuated, such that the momentum flux tendency of the easterly components of the FGW spectrum becomes small everywhere above 30 km. This, together with the decrease in the flux associated with westerly components of the spectrum noted above, produces a positive tendency in the net flux (green lines). Above 50 km, where the waves dissipate, this in turn gives rise to a negative tendency on FGWD, consistent with Fig. 6f.

5. Summary

WACCM4, the latest version of the Whole Atmosphere Community Climate Model, coupled to a fully interactive ocean model, has been used to investigate trends in the Brewer–Dobson circulation due to climate change. The model reproduces both branches of the climatological circulation as previously documented in ERA-Interim (e.g., Seviour et al. 2012) and other modeling studies (e.g., Birner and Bönisch 2011). The shallow branch, in the lowermost stratosphere, has upwelling in the tropics and downwelling in the subtropics; the deep branch, whose maximum upwelling is located in the upper stratosphere (~48 km) is global in extent, with downwelling extending to polar latitudes and down into the polar troposphere (not shown). A seasonal analysis shows the shift of the tropical pipe into the summer hemisphere in both branches of the circulation; this shift is more noticeable in the deep branch.

The application of the DCP between the turn-around latitudes of the climatological tropical upwelling revealed that different waves are involved in driving the climatology of the shallow and the deep branches. The climatological upwelling of the shallow branch is mainly driven by resolved waves (70%) with a nonnegligible contribution of OGWs (30%) at 70 hPa (~18.6 km), in agreement with Butchart et al.’s (2011) multimodel and reanalysis study. Resolved waves are also dominant in forcing the upwelling of the deep branch, followed by FGWs, which contribute negatively and consequently slow down the climatological circulation at these altitudes. Contributions from OGWs and CGWs are less than 10% to the deep branch upwelling.

Simulations using climate change scenarios show an acceleration of the BDC that increases with increasing loading of GHGs. Significant trends are found in both branches of the BDC in RCP4.5 and RCP8.5 scenarios (larger in RCP8.5), but the trends are much weaker and mostly not significant in the RCP2.6 scenario. Trends in tropical upwelling are confined to lower latitudes than the climatological upwelling in the lower and upper stratosphere, where the maximum upwelling in both branches occurs, which indicates a narrowing of the tropical upwelling region in the lower stratosphere and in the stratopause region. In addition, the maximum changes in the upwelling of the upper stratosphere appear at higher altitudes than the climatology, indicating an upward extension of the circulation in these climate change scenarios. In WACCM4, the largest changes in the tropical upwelling occur at 100 hPa and at 1 hPa. In the upper stratosphere, the tropical upwelling increases up to 50% in the RCP8.5 scenario compared to climatology over the period covered by this simulation.

To identify the waves involved in the acceleration of the mean meridional circulation, we have applied the downward control principle, choosing the climatological turn-around latitudes of the historical period as the limits of integration for calculating tropical and extratropical mean tendencies. The trend in the shallow branch upwelling is mainly due to waves that are explicitly resolved in the model (70%) while parameterized orographic gravity waves account for 30% in the lower stratosphere (at 70 hPa). This agrees with results from other chemistry–climate models such as the Atmospheric Model with Transport and Chemistry (AMTRAC) (Li et al. 2008), the Canadian Middle Atmosphere Model (CMAM) (McLandress and Shepherd 2009), and the previous version of WACCM, WACCM3.5 (GR08). In addition, Calvo and Garcia (2009) showed that the main driver of trends in the tropical upwelling of the lowermost stratosphere in WACCM3.5 was the EP flux divergence due to large-scale tropical Rossby waves. We note, nevertheless, that models are by no means unanimous in their attribution of trends in the shallow branch; in a minority of the models included in Fig. 4.11 of Eyring et al. (2010), the trend is due to changes in parameterized wave driving. In addition, the absolute magnitude of the trend, regardless of its source, varies considerably among models.

In the upper stratosphere, where frontal and convective gravity waves start to break, the contribution to the trend in the tropical upwelling is more heterogeneous than in the lower stratosphere. At these altitudes, frontal gravity waves play an important role in determining the upwelling trend, and their contribution becomes larger as the GHG burden increases. Trends in calculated from the DCP revealed that more than 50% of the tropical upwelling trend near the stratopause in simulation RCP8.5 is due to changes in frontal gravity wave drag while resolved wave drag accounts for 30%. The same methodology applied to the extratropics confirmed that frontal gravity waves are the main contributor to extratropical downwelling trends. The dominant role of frontal gravity waves in driving the mean meridional circulation trends in the upper stratosphere is related to changes in the background zonal-mean zonal winds. These changes in the winds modify the spectrum of frontal gravity waves such that the net momentum flux becomes more positive above the region where the subtropical jet intensifies (~15 km), such that the flux divergence in the upper stratosphere becomes more negative, which enhances tropical upwelling.

Finally, it should be borne in mind that the results presented here must be presumed to be model dependent. This is particularly the case in the upper stratosphere, where gravity waves are a major driver of the trends, since gravity waves are parameterized in general circulation models and these parameterizations differ among the models. Even in the lower stratosphere, Fig. 4.11 of Eyring et al. (2010) shows that, while all models produce positive tropical mass flux trends at 70 hPa, the magnitude of the trends and the relative contribution thereto by different waves varies from model to model. Nevertheless, insofar as the trends computed here depend on filtering of the (parameterized) upward-propagating gravity wave spectra by the evolving zonal-mean wind distribution, the sign and vertical distribution of the trends should be robust.

Acknowledgments

The authors are grateful to Jadwiga H. Richter for the helpful discussion and comments on gravity waves. F. M. Palmeiro and N. Calvo were supported by the Spanish Ministry of Economy and Competitiveness trough project CGL12012-34221, “Mechanisms and Variability of the Troposphere-Stratosphere Coupling.” This work was also partially supported by the European Project 603557 STRATOCLIM under program FP7-ENV-2013-two-stage. The Community Earth System Model (CESM) is supported by NSF and the Office of Science of the U.S. Department of Energy. Computing resources were provided by NCAR’s Climate Simulation Laboratory, sponsored by NSF and other agencies. This research was enabled by the computational and storage resources of NCAR’s Computational and Information Systems Laboratory (CISL).

REFERENCES

REFERENCES
Andrews
,
D. G.
,
J. R.
Holton
, and
C. B.
Leovy
,
1987
: Middle Atmosphere Dynamics. International Geophysics Series, Vol. 40, Academic Press, 489 pp.
Birner
,
T.
, and
H.
Bönisch
,
2011
:
Residual circulation trajectories and transit times into the extratropical lowermost stratosphere
.
Atmos. Chem. Phys.
,
11
,
817
827
, doi:.
Braesicke
,
P.
, and
J. A.
Pyle
,
2004
:
Sensitivity of dynamics and ozone to different representations of SSTs in the Unified Model
.
Quart. J. Roy. Meteor. Soc.
,
130
,
2033
2045
, doi:.
Butchart
,
N.
, and Coauthors
,
2006
:
Simulations of anthropogenic change in the strength of the Brewer–Dobson circulation
.
Climate Dyn.
,
27
,
727
741
, doi:.
Butchart
,
N.
, and Coauthors
,
2010
:
Chemistry–climate model simulations of twenty-first century stratospheric climate and circulation changes
.
J. Climate
,
23
,
5349
5374
, doi:.
Butchart
,
N.
, and Coauthors
,
2011
: Multimodel climate and variability of the stratosphere. J. Geophys. Res.,116, D05102, doi:.
Calvo
,
N.
, and
R. R.
Garcia
,
2009
:
Wave forcing of the tropical upwelling in the lower stratosphere under increasing concentrations of greenhouse gases
.
J. Atmos. Sci.
,
66
,
3184
3196
, doi:.
Danabasoglu
,
G.
,
S. C.
Bates
,
B. P.
Briegleb
,
S. R.
Jayne
,
M.
Jochum
,
W. G.
Large
,
S.
Peacock
, and
S. G.
Yeager
,
2012
:
The CCSM4 ocean component
.
J. Climate
,
25
,
1361
1389
, doi:.
Deckert
,
R.
, and
M.
Dameris
,
2008
:
Higher tropical SSTs strengthen the tropical upwelling via deep convection
.
Geophys. Res. Lett.
,
35
,
L10813
, doi:.
Eyring
,
V.
,
T. G.
Shepherd
, and
D. W.
Waugh
, Eds.,
2010
: SPARC report on the evaluation of chemistry climate models. SPARC Rep. 5, WCRP-132, WMO/TD-1526, 305 pp. [Available online at http://www.sparc-climate.org/publications/sparc-reports/sparc-report-no5/.]
Garcia
,
R. R.
, and
W. J.
Randel
,
2008
:
Acceleration of the Brewer–Dobson circulation due to increases in greenhouse gases
.
J. Atmos. Sci.
,
65
,
2731
2739
, doi:.
Garny
,
H.
,
M.
Dameris
,
W.
Randel
,
G. E.
Bodeker
, and
R.
Deckert
,
2011
:
Dynamically forced increase of tropical upwelling in the lower stratosphere
.
J. Atmos. Sci.
,
68
,
1214
1233
, doi:.
Hardiman
,
S.
,
N.
Butchart
, and
N.
Calvo
,
2014
:
The morphology of the Brewer–Dobson circulation and its response to climate change in CMIP5 simulations
.
Quart. J. Roy. Meteor. Soc.,
doi:, in press.
Haynes
,
P. H.
,
C. J.
Marks
,
M. E.
McIntyre
,
T. G.
Shepherd
, and
K. P.
Shine
,
1991
: On the “downward control” of extratropical diabatic circulations by eddy-induced mean zonal forces. J. Atmos. Sci., 48, 651–680, doi:.
Holland
,
M. M.
,
D. A.
Bailey
,
B. P.
Briegleb
,
B.
Light
, and
E.
Hunke
,
2012
:
Improved sea ice shortwave radiation physics in CCSM4: The impact of melt ponds and aerosols on Arctic sea ice
.
J. Climate
,
25
,
1413
1430
, doi:.
Holton
,
J. R.
,
1990
:
On the global exchange of mass between the stratosphere and troposphere
.
J. Atmos. Sci.
,
47
,
392
395
, doi:.
Holton
,
J. R.
,
P. H.
Haynes
,
M. E.
McIntyre
,
A. R.
Douglass
,
R. B.
Rood
, and
L.
Pfister
,
1995
:
Stratosphere–troposphere exchange
.
Rev. Geophys.
,
33
,
403
439
, doi:.
Kosaka
,
Y.
, and
S.-P.
Xie
,
2013
:
Recent global-warming hiatus tied to equatorial Pacific surface cooling
.
Nature
,
501
,
403
407
, doi:.
Lean
,
J.
,
G.
Rottman
,
J.
Harder
, and
G.
Kopp
,
2005
:
SORCE contributions to new understanding of global change and solar variability
.
Sol. Phys.
,
230
,
27
53
, doi:.
Li
,
F.
,
J.
Austin
, and
J.
Wilson
,
2008
:
The strength of the Brewer–Dobson circulation in a changing climate: Coupled chemistry–climate model simulations
.
J. Climate
,
21
,
40
–57, doi:.
Lin
,
P.
, and
Q.
Fu
,
2013
:
Changes in various branches of the Brewer–Dobson circulation from an ensemble of chemistry climate models
.
J. Geophys. Res. Atmos.
,
118
,
73
84
, doi:.
Marsh
,
D.
,
M. J.
Mills
,
D. E.
Kinnison
,
J.-F.
Lamarque
,
N.
Calvo
, and
L. M.
Polvani
,
2013
:
Climate change from 1850 to 2005 simulated in CESM1(WACCM)
.
J. Climate
, 26, 7372–7391, doi:.
Matthes
,
K.
,
D. R.
Marsh
,
R. R.
Garcia
,
D. E.
Kinnison
,
F.
Sassi
, and
S.
Walters
,
2010
: Role of the QBO in modulating the influence of the 11 year solar cycle on the atmosphere using constant forcings. J. Geophys. Res.,115, D18110, doi:.
McLandress
,
C.
, and
T. G.
Shepherd
,
2009
:
Simulated anthropogenic changes in the Brewer–Dobson circulation, including its extension to high latitudes
.
J. Climate
,
22
,
1516
–1540, doi:.
Olsen
,
M. A.
,
M. R.
Schoeberl
, and
J. E.
Nielsen
,
2007
:
Response of stratospheric circulation and stratosphere–troposphere exchange to changing sea surface temperatures
.
J. Geophys. Res.
,
112
,
D16104
, doi:.
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
, doi:.
Rosenlof
,
K. H.
,
1995
:
Seasonal cycle of the residual mean meridional circulation in the stratosphere
.
J. Geophys. Res.
,
100
,
5173
5191
, doi:.
Seviour
,
W. J. M.
,
N.
Butchart
, and
S. C.
Hardiman
,
2012
:
The Brewer–Dobson circulation inferred from ERA-Interim
.
Quart. J. Roy. Meteor. Soc.
,
138
,
878
888
, doi:.
Shepherd
,
T. G.
, and
C.
McLandress
,
2011
:
A robust mechanism for strengthening of the Brewer–Dobson circulation in response to climate change: Critical-layer control of subtropical wave breaking
.
J. Atmos. Sci.
,
68
, 784–797, doi:.
Tilmes
,
S.
,
R. R.
Garcia
,
D. E.
Kinnison
,
A.
Gettelman
, and
P. J.
Rasch
,
2009
: Impact of geoengineered aerosols on the troposphere and stratosphere. J. Geophys. Res.,114, D12305, doi:.

Footnotes

*

The National Center for Atmospheric Research is sponsored by the National Science Foundation.