The Brewer‐Dobson Circulation During the Last Glacial Maximum

The Brewer‐Dobson circulation during the Last Glacial Maximum (LGM) is investigated in simulations using the Whole Atmosphere Community Climate Model version 6. We examine vertical mass fluxes, age of stratospheric air, and the transformed Eulerian mean stream function and find that the modeled annual‐mean Brewer‐Dobson circulation during the LGM is almost everywhere slower than that in the modern climate (with or without anthropogenic ozone depleting substances). Compared to the modern climate, the annual‐mean tropical upwelling in the LGM is 11.3–16.9%, 11.2–15.8%, and 4.4–10.2% weaker, respectively, at 100, 70, and 30 hPa. Simulated decreases in annual‐mean mass fluxes at 70 and 100 hPa are caused by a weaker parameterized orographic gravity wave drag and resolved wave drag, respectively.


Introduction
The Brewer-Dobson circulation (BDC) is the global overturning mass circulation in which tropospheric air enters the stratosphere across the tropical tropopause, moves upward and poleward, and descends into the extratropical troposphere (e.g., Butchart, 2014;Haynes et al., 1991;Holton et al., 1995). The BDC varies with season, with a stronger cell occurring in the winter hemisphere where stratospheric winds are westerly (e.g., Plumb, 2002;Rosenlof, 1995). Because of considerably more longitudinal asymmetry in topography and land and ocean in the Northern Hemisphere (NH) than in the Southern Hemisphere (SH), planetary wave drag during winter is much stronger in the NH, leading to a stronger BDC.
The stratosphere during the LGM is expected to be considerably different from today because of reduced CO 2 , CH 4 , and N 2 O, the presence of large ice sheets in the NH, and increased latitudinal sea surface temperature (SST) gradients. Rind et al. (2001) used the GISS Climate Middle Atmosphere Model (8°× 10°a nd 23 vertical layers) to investigate how the stratosphere may have changed during the LGM. They suggested a warmer LGM stratosphere, due primarily to reduced CO 2 , and an increase in the BDC in the lower stratosphere but a decrease above. Using the GISS model with a higher resolution (4°× 5°and 53 layers), Rind et al. (2009) again found an intensified LGM BDC in the lower stratosphere and increased stratosphere-troposphere exchange (STE) at higher latitudes, but a decreased circulation for most of the middle atmosphere. Noting intensified BDC and STE associated with increased CO 2 (e.g., Butchart, 2014;Hegglin & Shepherd, 2009;Shepherd, 2008), Rind et al. (2009) argued that the LGM SSTs enhanced latitudinal temperature gradients, thus forcing a similar STE response. Here we investigate the BDC during the LGM using the Whole Atmosphere Community Climate Model version 6 (WACCM6). We find that the annual-mean BDC during the LGM is almost everywhere weaker than that in the modern climate.
The BDC controls the spatial distribution of stratospheric O 3 as well as stratosphere-troposphere O 3 transport, which have important implications for LGM tropospheric chemistry through the impact on ultraviolet radiation as well as a direct contribution to tropospheric O 3 (e.g., Geng et al., 2017;Murray et al., 2014). Quantifying the BDC in the LGM is important for interpretation of ice-core proxies of atmospheric chemistry (Geng et al., 2017;Yeung et al., 2019) and biological productivity (Luz et al., 1999). A dynamically consistent O 3 field in the LGM can also be expected to be important to the simulation of the LGM climate (e.g., Chiodo & Polvani, 2017;Nowack et al., 2015;Szopa et al., 2019). The present study should also broaden our perspective regarding the stratospheric response to anthropogenically induced climate changes in the future.

Model Experiments
We used the WACCM6, which is the high-top atmospheric component of the Community Earth System Model version 2, described by Gettelman et al. (2019). The model was run with specified ocean and ice components with a horizontal resolution of 0.9°longitude × 1.25°latitude. WACCM6 has 70 levels up to 140 km and includes the physics and chemistry of the troposphere, stratosphere, and mesosphere. The WACCM model suite has a realistic evolution of the SH springtime ozone hole over the latter half of the twentieth century  and a realistic frequency of stratospheric sudden warming (Gettelman et al., 2019). Further, as a result of improvements in the model climatology and the gravity wave parameterization, WACCM6 has an internally generated rather than prescribed QBO in the lower stratosphere (Garcia & Richter, 2019). We couple the WACCM6 model to the CLM4.0 land model and make use of the CLM4.0 LGM lower boundary conditions from the CESM Paleo Working Group (e.g., Brady et al., 2013).  (Rayner et al., 2003) with a seasonal evolution. For LGM PMIP3 , we used the PMIP3 multimodel mean SST/sea ice for all models that had both LGM and historical runs (CCSM4, CNRM-CM5, GISS-E2-R, IPSL-CM5A-LR, MIROC-ESM, MPI-ESM-P, and MRI-CGCM3) (Braconnot et al., 2012). The SSTs/sea ice distributions for the LGM were obtained from the model differences between LGM and current climate simulations plus current observed SSTs/sea ice. For LGM PROXY , the SSTs were from the MARGO proxy data (Kucera et al., 2005), while sea ice was derived as for LGM PMIP3 . The ice sheet topography from Abe-Ouchi et al. (2015) was used in the LGM simulations. The results from LGM PMIP3 and LGM PROXY inform the sensitivity to the SSTs used for the LGM.
Before investigating the BDC during the LGM, we first examine WACCM6's performance by comparing the WACCM6 MC simulation with the ERA-I reanalysis and quantify the ODS effects on the BDC (section 3).

MC Simulation Versus ERA-I Reanalysis and ODS Effect
Supporting information Table S1 shows the annual-mean air mass fluxes at different altitudes and by different BDC branches from the MC simulation compared with ERA-I reanalysis. The BDC was decomposed into LGM PMIP3 and MC. (Right) relative differences in mass fluxes transported by transition (100-70 hPa) and stratospheric shallow (70-30 hPa) and deep (above 30 hPa) branches of the Brewer-Dobson circulation. Black, red, and blue colors indicate results associated with upwelling in the tropics and downwelling in the Northern Hemisphere and Southern Hemisphere, respectively. The annual-mean and seasonal results are shown in row (a) and rows (b-e), respectively. The relative differences that are statistically significant at 95% confidence level are denoted with the slant lines in the bar for the middle and left panels.

Geophysical Research Letters
the transition (100-70 hPa) and stratospheric shallow (70-30hPa) and deep (30 hPa-) branches following Lin and Fu (2013). The WACCM6 MC compares well with ERA-I and shows better agreement with the reanalysis than the average of 12 CCMs (see Table 1 of Lin & Fu, 2013). The MC simulation also shows excellent agreement with ERA-I for the spatial distribution of annual-mean residual vertical velocities at 70 hPa ( Figure S1) and annual-and seasonal-mean transformed Eulerian mean (TEM) stream functions ( Figure S2). The annual-mean age of air for the MC ( Figure S3) also agrees well with previous estimates (e.g., Garcia et al., 2011;Li et al., 2012).
Compared with the MC (Table S2), the MC no_ODS annual-mean upwelling mass fluxes decrease by 1.8% and 3.6%, respectively, at 70 and 30 hPa (Table S3), which are both statistically significant. The mass flux transported by the SH deep branch in DJF (December, January, and February) is only 0.44 × 10 9 kg/s in MC (Table S2) but decreases by 39.8% in MC no_ODS (Table S3), which is largely responsible for simulated significant changes in the annual-mean upwelling mass fluxes at 70 and 30 hPa. This result agrees with Lin and Fu (2013) that the response of the BDC to the ozone depletion is primarily in the southern deep branch during the austral summer. Further, significant changes are also seen in SH transition and stratospheric shallow branches in various seasons (Table S3).

BDC During the LGM Versus That in the Modern Climate
The far left panel in Figure 1 shows annual-and seasonal-mean mass flux profiles from the MC and LGM PMIP3 for the upwelling in the tropics and its partitioning into downwelling in the NH and SH. While the mass fluxes show large dependences on season, hemisphere, and altitude, the fluxes from LGM PMIP3 are almost always smaller than those from MC. Shown in the middle panel are the corresponding relative differences in mass fluxes at 100, 70, and 30 hPa between LGM PMIP3 and MC: The annual-mean tropical upwelling from LGM PMIP3 is 16.9%, 15.8%, and 7.8% smaller than those in MC, respectively. Most of the relative differences shown in the middle panels of Figure 1 for the four seasons are also negative and statistically significant. Except for the SH in DJF, stratospheric ozone depletion has a small impact (compare Figure 1 with Figure S4, in which MC has been replaced with MC no_ODS ). The annual mean tropical upwelling fluxes at 100, 70, and 30 hPa from LGM PMIP3 are 16.5%, 14.2%, and 4.4% smaller than those from MC no_ODS , respectively. The far right panel of Figure 1 shows the relative differences in mass fluxes transported by the BDC transition and stratospheric shallow and deep branches for LGM PMIP3 versus MC. The stratospheric ozone depletion effects (compare Figures 1 and S4) are again associated with the stratospheric deep cell in the SH in DJF. The annual-mean mass fluxes transported by the BDC transition and stratospheric shallow and deep branches from LGM PMIP3 are 18.4% (19.3%), 23.3% (23.2%), and 7.8% (4.4%) smaller than those from MC (MC no_ODS ), respectively, and these changes are all statistically significant.
The age of stratospheric air is defined as the time to transport an air parcel from the tropical tropopause to a given location in the stratosphere, which is an indicator of the strength of the BDC (Hall & Plumb, 1994;Waugh & Hall, 2002). The age of stratospheric air is derived in WACCM6 by using a synthetic LGM PMIP3 versus MC no_ODS , indicating increased age of air everywhere in the stratosphere, as would be expected based on the decrease in the annual-mean BDC during the LGM (top panels in Figures 1 and S4). The larger increases in (a) compared to (b) are caused by stratospheric ozone depletion. The results shown in Figure 2 are somewhat different from those in Rind et al. (2009) (see their Figure 3). Although Rind et al. (2001Rind et al. ( , 2009 also showed a decrease in the middle atmosphere residual circulation (e.g., over~50 hPa), they found an increase in the lower stratosphere circulation over the extratropics.
The Lagrangian-mean meridional mass transport by the BDC can be approximated by the TEM stream function (Andrews et al., 1987). Figure 3 shows the annual-and seasonal-mean TEM stream functions from the MC simulation (contours) and the differences between LGM PMIP3 and MC (colors). Positive (negative) values indicate clockwise (counterclockwise) turning. For the region with solid (dashed) contours, a negative (positive) change indicates a decrease of the circulation intensity. We see a deceleration of the annual-mean BDC in both the NH and SH cells (Figure 3a). The seasonal-mean values of the BDC strength also show an overall deceleration (Figures 3b and 3e) but with a few acceleration regions including parts of the NH (SH) cell in DJF (JJA) and parts of both NH and SH cells in MAM and SON. Similar results are obtained by comparing the LGM PMIP3 with the MC no_ODS simulations ( Figure S5). It is interesting to compare the TEM stream function changes in DJF and JJA (Figures 3b and 3d) to those obtained by Rind et al. (2001) (see the left panel in their Plate 6). In DJF, the BDC from Rind et al. (2001) is largely weakening in the middle stratosphere, consistent with Figure 3b, but it is strengthening in the lower stratosphere below~50 hPa, especially in the NH, which is different from Figure 3b. In JJA, the changes in the TEM stream function from Rind et al. (2001) show a largely strengthening BDC in the SH, and also a strengthening in the NH between 200 and 100 hPa, different from Figure 3d. The global-mean surface air temperatures are 15.0, 10.5, and 10.0°C from the WACCM6 MC, LGM PROXY , and LGM PMIP3 simulations, respectively, corresponding to global-mean SSTs of 18.4, 16.6, and 16.2°C. The global-mean surface air temperatures and SSTs from LGM PROXY are thus~0.4-0.5°C higher than LGM PMIP3 . The higher global-mean SSTs from LGM PROXY are due to higher tropical SSTs compared to LGM PMIP3 ; middle-latitude SSTs are lower in LGM PROXY ( Figure S6). This leads to a larger latitudinal SST gradient in LGM PROXY . The BDC from LGM PROXY is generally stronger than that from LGM PMIP3 , especially for the transition and stratospheric shallow branches (Table S4). However, the results from LGM PROXY relative to MC and MC no_ODS (Figures S7-S11) are very similar to results from LGM PMIP3 relative to MC and MC no_ODS (Figures 1-3, S4, and S5). The annual-mean tropical upwelling at 100, 70, and 30 hPa from LGM PROXY is 11.7 (11.3)%, 12.8 (11.2)%, and 10.2 (6.9) % smaller than those from MC (MC no_ODS ), respectively ( Figures SM7 and SM10).

Geophysical Research Letters
FU ET AL. Figure 4 shows the annual-mean air mass fluxes derived from the downward control principle versus those from the TEM residual vertical velocity at 70 and 100 hPa for tropical upwelling and NH and SH downwelling, for the four WACCM6 simulations. The mass fluxes from the two methods agree well. The decrease in the modeled mass fluxes at 70 hPa during the LGM is largely caused by weaker parameterized orographic gravity wave drag (Figure 4a). This is because of a downward shift of the NH subtropical jet and slow down of the SH jets during the LGM (Figure 5a), both of which cause the parameterized orographic gravity waves to break lower in the atmosphere (Figures 5b and 5d) (Butchart, 2014;Butchart et al., 2010;Li et al., 2008;Okamoto et al., 2011). This results in a less drag above 70 hPa, which leads to weaker tropical upwelling at that level. Note that following the downward control principle the residual vertical velocity at a given altitude is determined by the vertically integrated wave forcing above that altitude and changes in the wave drag in the vicinity of the turnaround latitudes are most relevant. Our result here is broadly consistent with Butchart et al. (2010), who found that the percentage contribution of the parameterized waves to driving trends in the 70 hPa upwelling caused by an increase in CO 2 is much larger than their percentage contribution to driving the upwelling itself. Similar to the parameterized gravity waves, a downward shift of wave dissipation is also found in the resolved waves (Figures 5c and 5e) due to the changes in the subtropical jets ( Figure 5a) (Cohen et al., 2014;Garcia & Randel, 2008;Shepherd & McLandress, 2011). Less wave drag from the resolved waves is largely responsible for the mass flux decreases at 100 hPa (Figures 4b, 5c, and 5e). The results in this paper are not sensitive to the subgrid-scale topography treated in the LGM simulations (e.g., the present day versus LGM subgrid-scale topography) (not shown). Rind et al. (2009) found a larger global poleward eddy heat flux in their LGM simulations than their present-day control run, which was responsible for the strengthening of their LGM BDC in the lower stratosphere. Based on our simulations (Figure 4), although the contribution of resolved waves to the mass fluxes at 70 hPa during the LGM also becomes slightly larger in the NH than the modern climate, it is not enough to compensate for the decrease caused by the orographic gravity wave drag changes. It is also interesting to note that in Rind et al. (2001), the annual-mean surface air temperature during LGM based on the proxy data is 4.3°C colder than in their control run, which is very similar to our 4.5°C from LGM PROXY versus MC. Rind et al. (2001Rind et al. ( , 2009 showed that their results are not sensitive to either SSTs used or wave drag treatments associated with the ice sheet. The differences between our results and Rind et al. (2001Rind et al. ( , 2009) highlight the need for multimodel simulations and comparison of the LGM stratospheric circulations.

Summary and Conclusions
This study examines the BDC in the LGM with simulations for both LGM and modern climate using WACCM6. The BDC simulation for the modern climate agrees well with ERA-I in terms of annual-mean air mass fluxes at different levels and by different BDC branches as well as the annual-and seasonal-mean TEM stream functions.
We show smaller annual-mean air mass fluxes during the LGM compared to the modern climate at all levels and branches. The annual-mean tropical upwelling in the LGM is 11.3-16.9%, 11.2-15.8%, and 4.4-10.2% smaller, respectively, at 100, 70, and 30 hPa than in the modern climate (with or without anthropogenic ozone depleting substances). The age of stratospheric air in the LGM increases everywhere as compared with that from the modern climate. The annual-mean TEM stream functions also indicate a slowdown of the BDC during the LGM. The annual-mean BDC during the LGM is therefore always slower than that in the modern climate regardless of altitude, latitude, and BDC vertical and hemispheric branches. This result is insensitive to the SSTs used in the LGM simulation and irrespective of the presence of anthropogenic ODS in the modern climate. The seasonal-mean BDC during the LGM also shows an overall slowdown of the circulation apart from a few regions.
A decreasing BDC in the LGM will have important implications for the spatial distribution of ozone in the stratosphere, as well as STE and surface ultraviolet radiation. The latter two changes can be expected to affect tropospheric oxidant abundances, with potential implications for the lifetime of trace gases such as methane. These impacts will be explored in future studies.