Sensitivity of simulated convection‐driven stratosphere‐troposphere exchange in WRF‐Chem to the choice of physical and chemical parameterization

Abstract

Tropopause‐penetrating convection is capable of rapidly transporting air from the lower troposphere to the upper troposphere and lower stratosphere (UTLS), where it can have important impacts on chemistry, the radiative budget, and climate. However, obtaining in situ measurements of convection and convective transport is difficult and such observations are historically rare. Modeling studies, on the other hand, offer the advantage of providing output related to the physical, dynamical, and chemical characteristics of storms and their environments at fine spatial and temporal scales. Since these characteristics of simulated convection depend on the chosen model design, we examine the sensitivity of simulated convective transport to the choice of physical (bulk microphysics or BMP and planetary boundary layer or PBL) and chemical parameterizations in the Weather Research and Forecasting model coupled with Chemistry (WRF‐Chem). In particular, we simulate multiple cases where in situ observations are available from the recent (2012) Deep Convective Clouds and Chemistry (DC3) experiment. Model output is evaluated using ground‐based radar observations of each storm and in situ trace gas observations from two aircraft operated during the DC3 experiment. Model results show measurable sensitivity of the physical characteristics of a storm and the transport of water vapor and additional trace gases into the UTLS to the choice of BMP. The physical characteristics of the storm and transport of insoluble trace gases are largely insensitive to the choice of PBL scheme and chemical mechanism, though several soluble trace gases (e.g., SO₂, CH₂O, and HNO₃) exhibit some measurable sensitivity.

Introduction

Extratropical convection, capable of penetrating the tropopause and reaching up to 20 km in altitude, plays an important role in the exchange of trace gases, such as ozone and water vapor, between the stratosphere and troposphere (stratosphere‐troposphere exchange or STE) [Holton et al., ]. The rapid transport of these gases [Mullendore et al., ], many of which typically have a short lifetime near the surface, can have a large impact on the chemistry of the atmosphere when transported to the upper troposphere [Barth et al., ]. Specifically, transported species, such as nonmethane hydrocarbons (NMHCs), peroxides, formaldehyde, and CH₃OH, can react to form HO₂ and RO₂ and produce ozone [e.g., Pickering et al., ; Barth et al., ] while lightning‐produced nitrogen oxides (NO_x) are also an important precursor in ozone production [Ridley et al., , ]. NMHCs, SO₂, and volatile organic compounds (VOCs) can react and create new aerosols when transported to the UT [e.g., Thornton et al., ]. Some halogen species can be transported to the stratosphere where they can affect ozone chemistry [Dvortsov et al., ]. Additionally, exchange of longer‐lived greenhouse gases such as ozone and water vapor within convection and from other processes can significantly impact the radiation budget [e.g., Forster and Shine, ]. However, despite these important impacts, the effects of tropopause‐penetrating convection on the transport and exchange of trace gases are not well understood due to limitations of existing observing systems and computational efficiencies that allowed atmospheric chemistry simulations to be run at horizontal resolutions finer than 4 km.

Many recent observational and modeling studies have focused on documenting the effects of extratropical tropopause‐penetrating convection on the composition of the UTLS and investigating the mechanisms responsible for exchange. Aircraft observations from individual events have been used to identify how tropopause‐penetrating convection redistributes trace gases in the atmosphere [e.g., Fischer et al., ; Hegglin et al., ; Ray et al., ; Hanisco et al., ; Anderson et al., ; Homeyer et al., ; Pan et al., ]. However, in situ observations of the convective overshoot (that extending well above the tropopause) are not possible with current research aircraft. Thus, modeling studies have been employed to provide new insights that could not be found through in situ observations [e.g., Stenchikov et al., ; Gray, ; Wang, ; Mullendore et al., ; Lane and Sharman, ; Bigelbach et al., ]. In particular, these studies have identified turbulent mixing from gravity wave and Kelvin‐Helmholtz instabilities (i.e., wave breaking) as mechanisms for irreversible exchange between the troposphere and stratosphere. Recent aircraft observations have revealed an additional and potentially unique transport mechanism: wrapping of stratospheric air around the perimeter of detraining anvil clouds and into the upper troposphere, which has been reproduced in a numerical model [Pan et al., ]. The dynamical process(es) responsible for this anvil wrapping mechanism are not yet known.

To better understand the mechanisms responsible for transport and mixing and the significance of convection‐driven transport to UTLS composition, new numerical simulations of extratropical tropopause‐penetrating convection are needed. However, since the evolution, vertical extent, and intensity of convection are sensitive to the model design, it is important to determine the choices that best reproduce the physical and chemical transport characteristics of observed storms. Model sensitivity tests for the choice of horizontal and vertical grid resolution have recently been completed and show that the depth of overshooting and cross‐tropopause transport increase with finer horizontal grid spacing but decrease with finer vertical grid spacing [Homeyer, ]. This paper will present results of model sensitivity to the choice of physical and chemical parameterization.

Model Design, Data, and Methods

Model Description and Initialization

Version 3.7.1 of the Weather Research and Forecasting model [Skamarock et al., ] coupled with Chemistry [Grell et al., ; Fast et al., ] is used in this study. Simulations are run with one‐way nesting from a parent domain with a horizontal resolution of 10 km to a nested domain with 2 km resolution. The vertical grid consists of 101 levels with a nominal grid spacing of ~250 m in the free troposphere and a model top of 30 hPa (~24 km). A 5 km damping layer is employed to prevent reflection of spurious waves off the model top. Meteorological initial and boundary conditions are provided every 6 h from the ERA‐Interim reanalysis, which are available with a horizontal resolution of ~80 km and a vertical resolution ranging from 650 to 1000 m in the extratropical UTLS [Dee et al., ]. Chemical initial and boundary conditions are defined using output from the Model of Ozone and Related chemical Tracers, version 4 (MOZART‐4) [Emmons et al., ].

Three bulk microphysics parameterizations (BMP), planetary boundary layer (PBL) parameterizations, and chemical mechanisms are tested in this study. The full list of simulations conducted is provided in Table . The chosen microphysics parameterizations are the Morrison 2‐moment (MOR) [Morrison et al., ], the Milbrandt and Yau 2‐moment (MY) [Milbrandt and Yau, ], and the NSSL 2‐moment (NSSL) [Mansell et al., ] schemes. The chosen PBL parameterizations are the Yonsei University (YSU) [Hong et al., ], the quasi‐normal scale elimination (QNSE) [Sukoriansky et al., ], and the asymmetric convective model, version 2 (ACM2) [Pleim, ] schemes. The chosen chemical mechanisms are the Regional Atmospheric Chemistry Model (RACM‐ESRL) [Stockwell et al., ; Ahmadov et al., ] coupled with the Modal Aerosol Dynamics Model for Europe/Secondary Organic Aerosol Model (MADE/SORGAM) [Ackermann et al., ; Schell et al., ], the Carbon Bond Mechanism, version Z, (CBMZ) [Zaveri and Peters, ] coupled with Model of Simulating Aerosol Interactions and Chemistry 4‐bin aerosol model (MOSAIC) [Zaveri et al., ], and the Model of Ozone and Related chemical Tracers chemical mechanism (MOZART) [Emmons et al., ] coupled with the Georgia Institute of Technology—Goddard Global Ozone Chemistry Aerosol Radiation and Transport aerosol model (GOCART) [Chin et al., ; Ginoux et al., ]. While testing the model sensitivity to each type of parameterization, the remaining two were held constant, resulting in seven independent simulations of each targeted convective event (described in more detail below). In particular, the schemes held fixed in each sensitivity test are the NSSL microphysics parameterization, YSU PBL parameterization, and RACM‐ESRL chemical mechanism.

Description of Simulations Performed

Test Group	BMP	PBL	CHEM
Control run	NSSL	YSU	RACM‐ESRL
BMP‐1	MOR	YSU	RACM‐ESRL
BMP‐2	MY	YSU	RACM‐ESRL
PBL‐1	NSSL	QNSE	RACM‐ESRL
PBL‐2	NSSL	ACM2	RACM‐ESRL
CHEM‐1	NSSL	YSU	MOZCART
CHEM‐2	NSSL	YSU	CBMZ

Additional model design choices held constant in all simulations are as follows: Smagorinsky first‐order closure was used for horizontal subgrid‐scale mixing. The RRTMG scheme was used for both short‐wave and long‐wave radiation [Iacono et al., ]. Anthropogenic emissions were generated using the 2011 National Emissions Inventory. Biogenic emissions were calculated online with the Model of Emissions of Gases and Aerosols from Nature (MEGAN V2.04) [Guenther et al., ]. Photolysis rates were calculated using the Fast‐J scheme [Wild et al., ; Fast et al., ].

The reflectivity calculation for the simulations was done using the Center for Analysis and Prediction of Storms (CAPS) Polarimetric Radar data Simulator (CAPS‐PRS) [e.g., Jung et al., ; Dawson et al., ]. The CAPS‐PRS is capable of computing polarimetric radar variables from high‐resolution numerical simulations using microphysics schemes with one, two, and three moments. Polarimetric variables are expressed as a function of hydrometeor mixing ratio as well as their drop size distribution and densities. The melting layer is accounted for using a continuous melting process for the entire spectrum of varying density and dielectric constants [Jung et al., ]. However, the simulator may overpredict simulated reflectivity in the stratiform rain region due to the assumption of a fixed DSD intercept parameter for hail [Jung et al., ]. Reflectivity calculations were done using a T‐matrix method, which differentiates between Rayleigh and Mie scattering where appropriate and is more accurate than the online WRF radar reflectivity calculation, which only uses Rayleigh scattering.

Descriptions of the Model Parameterizations

Tables  list some distinguishing features of each parameterization tested in this study. Justification for choosing each scheme is briefly elaborated upon here.

Predicted Hydrometeor Moments and Other Unique Details of the Bulk Microphysics Parameterizations (BMPs) Employed in This Study

BMP	Predicted Mass Mixing Ratios	Predicted Number Concentrations	Other
MOR	Qc, Qr, Qs, Qi, Qg/h	Nc, Nr, Ns, Ni, Ng/h	Prognostic aerosols turned on; hail option turned on
NSSL	Qc, Qr, Qs, Qi, Qg, Qh	Nc, Nr, Ns, Ni, Ng, Nh	Prognostic aerosols turned on; predicts graupel and hail density; assigns shape parameter based on graupel/hail density
MY	Qc, Qr, Qs, Qi, Qg, Qh	Nc, Nr, Ns, Ni, Ng, Nh	Allows shape parameter to vary in particle size distribution

Bulk Microphysics

Storms that penetrate the tropopause are often very intense convective events that produce large hail, strong winds, and sometimes tornadoes at the Earth's surface. These storms also have unique physical features, such as strong updrafts, overshooting tops, and sometimes above‐anvil cirrus plumes, which are important for STE of water vapor. Since the physical characteristics of these storms depend on microphysical processes, model representations can vary considerably with the choice of BMP. In addition, since each BMP assumes certain hydrometeor types and sedimentation rates, we expect the choice of BMP to impact the simulated hydration of the UTLS. Based on these anticipated sensitivities, three BMPs of second order or higher were selected (Table ) instead of single‐moment schemes, since double‐moment schemes are generally more successful at reproducing observed convective storms [Igel et al., ].

The MOR, MY, and NSSL 2‐moment schemes were chosen. These schemes differ in the number of predicted hydrometeor classes, assumptions about the gamma function in the size distribution, and the ability to interact with chemistry modules. MY and NSSL each predict the number concentration and mass mixing ratio for six hydrometeor classes: cloud droplets, rain, snow, cloud ice, graupel, and hail. MOR only predicts five hydrometeor classes: cloud droplets, rain, snow, cloud ice, and graupel or hail. In this study, hail is predicted for MOR. In the gamma function, MY allows the shape parameter to vary as a function of mean‐mass diameter, while the other two BMPs do not. NSSL, however, allows the fall speeds of graupel and hail to vary based on the density and through adjustments to the drag coefficient. Lastly, MOR and NSSL allow prognostic aerosols, which allow WRF‐Chem to compute the increase in number concentration of droplets with the BMP then handling condensation. However, while prognostic aerosol effects were turned on, other aerosol effects (aerosol‐radiation) were turned off.

Planetary Boundary Layer

The PBL scheme in a numerical model such as WRF‐Chem determines the subgrid‐scale turbulence and vertical mixing, as well as the vertical thermodynamic and kinematic profiles of the lowest model levels. Thus, the choice of PBL scheme will impact the simulated PBL height, the location and timing of convective initiation, and the source of the air ingested by the storm and transported to the UTLS. Current PBL schemes available in WRF‐Chem can be classified into three different groups based on their treatment of subgrid‐scale mixing: local, nonlocal, and hybrid (i.e., both local and nonlocal). The basic difference between local and nonlocal mixing involves the depth over which vertical levels influence a variable at a given point. That is, in a local mixing scheme, a variable is only influenced by variables at directly adjacent vertical levels. In a nonlocal mixing scheme, a variable may be influenced by a number of vertical levels. For a hybrid scheme, the stability of the boundary layer determines whether local or nonlocal mixing is used (i.e., for neutral or stable conditions, it uses local closure and turns off nonlocal transport). To evaluate PBL schemes here, we chose one scheme from each category (Table ): QNSE (local), YSU (nonlocal), and ACM2 (hybrid).

Characteristics of the Planetary Boundary Layer (PBL) Parameterizations Employed in This Study

	Mixing	Advantages	Disadvantages
YSU	Nonlocal (meteorology) and local (chemistry)	Simulates deeper vertical mixing in buoyancy‐driven PBL	Too deep PBL for springtime convection;
QNSE	Local	Intended to account for wave phenomena within stable PBL	Too cool, moist, and shallow PBL
ACM2	Hybrid	Potential temperature and velocity depicted with greater accuracy; same mixing for chemical species; and default BL model in CMAQ	Too deep PBL

Adapted from Cohen et al. [].

Chemical Mechanism

In WRF‐Chem, numerous options are available for the chemical mechanism. Among the 30 different options are seven gas‐phase chemical mechanisms with different couplings to aerosol schemes, optional incorporation of aqueous chemistry, adaptions to use the Kinetic Preprocessor library and the Rosenbrock solver [Sandu and Sander, ], or other deviations from the default procedure of each mechanism. Since the chemical effects of convection‐driven STE can be represented differently depending on the chemical mechanism used, three commonly used gas‐phase mechanisms were chosen (Table ), each coupled to aerosol modules of varying complexity.

Characteristics of the Chemical Mechanisms Employed in This Study

Chemistry	RACM‐ESRL	CBMZ	MOZART
Aerosol scheme	MADE/SORGAM/VBS (Modal—3 modes)	MOSAIC (Bin—4 bins)	GOCART (Bulk)
Aerosol/cloud interactions	Yes	Yes	No
SOA formation	Ahmadov et al. []	Zaveri et al. []; Hodzic and Jimenez []	None
Wet scavenging	Easter et al. []	Easter et al. []	Neu and Prather []
No. reactions	214 (gas‐phase)	132	157 (gas‐phase)
No. species	77	52	85

The representation of aerosol chemistry varies from a simple representation in MOZCART to more complex approaches in CBMZ‐MOSAIC and RACM‐ESRL with MADE/SORGAM. In particular, RACM‐ESRL is coupled to the MADE/SORGAM/VBS aerosol module that represents secondary organic formation (SOA) and thus is designed to give more representative concentrations of particulate matter at diameters 2.5 μm and less (PM_2.5) than the other two chemical mechanisms presented here. SOA formation in CBMZ‐MOSAIC follows a simplified approach [e.g., Zaveri et al., ; Hodzic and Jimenez, ] and MOZCART does not include SOA formation. The SOA formation used in CBMZ‐MOSAIC follows an empirical parameterization based on the ratio between observed SOA concentrations to excess CO and the photochemical age of the air mass. In this formulation, organic mass is emitted as lumped SOA precursor surrogate in proportion to anthropogenic or biomass burned CO emissions following the observed SOA and CO ratio in aged air. This surrogate then reacts with OH to form a single nonvolatile species that condenses to form SOA [Hodzic and Jimenez, ]. The SOA parameterization in RACM‐ESRL follows the method described in Ahmadov et al. []. Briefly, the parameterization is based on a four bin volatility basis set. VOCs are oxidized by the hydroxyl radical, ozone, and nitrate radical into anthropogenic and biogenic compounds. Organic mass in each bin is produced for both high and low NO_x regimes and partitioned into aerosol and gas phase. Couplings to less complex aerosol modules (i.e., fewer bins or modes) were chosen over more detailed treatments since the main goal of this study is to assess the model sensitivity for bulk transport of trace gases and water. The wet scavenging scheme employed in RACM‐ESRL and CBMZ follows the Easter et al. [] approach while in MOZCART follows the Neu and Prather [] approach for gases and Easter et al. [] for aerosols. Both methods treat wet deposition by grid‐resolved precipitation, scavenging of cloud‐phase aerosols and gases by collection and freezing, interstitial‐phase aerosols by impaction, and gas‐phase gases by mass transfer and reaction. The main difference between the two approaches is the treatment of HNO₃. In the Neu and Prather [] approach, HNO₃ is partitioned into cloud ice as a function of temperature based on a burial model. As noted in section 2.2.1, prognostic aerosols were turned on for these simulations. However, aerosol‐cloud interactions only worked with the RACM‐ESRL and CBMZ scheme, and not MOZCART. Additionally, aerosol effects on the radiation scheme were turned off.

Observations for Model Evaluation

To conduct a robust evaluation of the selected model parameterizations, three observed cases from the Deep Convective Clouds and Chemistry (DC3) field campaign [Barth et al., ] were simulated: 19–20 May 2012, 29–30 May 2012, and 1–2 June 2012. Figure  shows the DC3 flight paths for each case and the boundary of the WRF nested domain used for analysis. For brevity, we only include results from the 19–20 May case in this paper, but the results for the remaining cases are comparable (see supporting information Figures S1–S16).

Trace gas measurements from two DC3 aircraft are used to evaluate the model simulated trace gas distributions in the troposphere and lower stratosphere: the National Science Foundation‐National Center for Atmospheric Research (NSF‐NCAR) Gulfstream V (GV) and the National Aeronautics and Space Administration (NASA) DC‐8. Measurements of ozone (O₃), carbon monoxide (CO), and water vapor (H₂O) from each aircraft were obtained at a rate of 1 Hz, which corresponds to a horizontal resolution of 100–200 m at aircraft cruise speed. Measurements of formaldehyde (CH₂O), sulfur dioxide (SO₂), and nitric acid (HNO₃) were obtained at a resolution of 1–2 s, 2–10 s, and 1–2 s, respectively. Cloud measurements from the NASA DC‐8 were incomplete during DC3 especially before 29 May when the cloud particle imager was added to the aircraft, such that a cloud indicator based on forward‐facing video camera and cloud particle imaging (when available) was created to facilitate analysis of cloudy and clear‐sky conditions during the entire campaign. Additional detail on this manually developed cloud indicator is available in Pollack et al. []. On the NSF‐NCAR GV, cloud particles were measured throughout the field campaign using the 2D‐C probe [National Center for Atmospheric Research, ]. More information about the instruments on board the GV and DC‐8 used in this study is given in Table .

Precision and Uncertainty of Measurements From Aircraft‐Based Instruments Used in Model Evaluation

Instrument	Variable	Precision and Uncertainty	Aircraft	Reference
NCAR	O3	2 ppbv ± 5%	GV	Ridley et al. []
NCAR vUV	CO	3 ppbv ± 3%	GV	Similar to Gerbig et al. []
VCSEL	H2O	<2% ± 5%	GV	Zondlo et al. []
GT‐CIMS	SO2	0.0119 ppb	GV	Kim et al. []
GT‐CIMS	HNO3	0.0396 ppb	GV	Huey []
CAMS	CH2O	0.02–0.03 ppb ± 0.02–0.03 ppb	GV	Fried et al. []
2D‐C	Cloud	N/A	GV	NCAR []
NOAA CL	O3	0.6 ppbv ± 5%	DC8	Davis et al. []; Dorsi et al. []
DACOM	CO	2 ppbv ± 2%	DC8	Ryerson et al. []
NASA DLH	H2O	1 ppmv ± 5%	DC8	Diskin et al. []
GT‐CIMS	SO2	0.002 ppb	DC8	Kim et al. []
GT‐CIMS	HNO3	0.1 ppb ± 50%	DC8	Huey []
DFGAS	CH2O	0.05 ppb ± 0.05 ppb	DC8	Fried et al. []
Cloud indicator	Cloud	‐‐‐	DC8	‐‐‐

The structure, intensity, and evolution of the simulated storms are evaluated in this study using three‐dimensional composites of individual radar volumes from the Next Generation Weather Radar (NEXRAD) program Weather Surveillance Radar—1988 Doppler (WSR‐88D) network [Crum and Alberty, ]. When operating in convective mode, data from individual radars are available every 4–7 min on native spherical grids with 14 elevation scans. These individual radar volumes were obtained from the National Centers from Environmental Information and merged into large‐area multiradar composites following the methods outlined in Homeyer [] and updated in Homeyer and Kumjian []. The composites have a temporal resolution of 5 min, a horizontal resolution of 0.02° longitude latitude (~2 km), and a vertical resolution of 1 km.

Background and Analysis of 19 May DC3 storm

On 19 May 2012 a deep convective line in central Oklahoma and Kansas initiated along a surface cold front. The storm was sampled by the NSF‐NCAR GV and NASA DC‐8 between 21:00 UTC on 19 May and 02:15 UTC on 20 May. A double tropopause was present at the time as a result of a poleward breaking Rossby wave, which transported air from the tropical upper troposphere to the extratropical lower stratosphere [e.g., Pan et al., ; Homeyer et al., ]. It has also been hypothesized that the double tropopause environment facilitated the large depth of overshooting, with storm top altitudes reaching up to 4 km above the unperturbed primary tropopause [Homeyer et al., ].

The analysis period was 22:30 UTC–03:00 UTC for the model simulations and 21:00 UTC–02:15 UTC for the observed storm, corresponding to comparable stages in the evolutions of the observed and modeled storms. The analysis area was a subset of the nested domain shown in Figure , which was subjectively determined in the model simulations to correspond to the convective system sampled by the DC3 aircraft. All figures shown were generated using this analysis area and time frame, unless otherwise noted.

Results and Discussion

Structure and Organization of the Simulated 19 and 20 May Storm

Figure  shows simulated column‐maximum radar reflectivity for each independent model run and the observed storm at comparable times.

Based on the radar reflectivity fields, it is apparent that the simulated storm is most sensitive to the choice of BMP and least sensitive to the choice of chemical mechanism. While the choice of BMP has little impact on the timing and location of the simulated storm (Figure S20), there are some clear differences in the horizontal scale and mode of convective organization. In particular, the NSSL BMP gives a more horizontally narrow storm than MOR or MY, which best agrees with the observed storm. This is likely due to the NSSL scheme's design, which allows the density of hail and graupel to vary and arguably leads to more reliable differential size sorting and sedimentation rates as described in Dawson et al. []. In the study performed by Dawson et al. [], idealized simulations of an observed supercell storm with a 3‐moment version of the NSSL bulk microphysics scheme were done while allowing and disallowing size sorting for hydrometeor species, considering several velocity‐diameter relationships for hail fall speed, and compared fixed and variable bulk densities that span the graupel‐to‐hail spectrum. Their simulations showed that the best performing simulations were the ones in which size sorting was allowed for rain and hail, and the bulk density and fall‐speed curve for hail were predicted. Thus, the differential size sorting of the NSSL BMP is likely the reason for the better representation of the simulated storm.

In contrast to the BMPs, the physical characteristics of the simulated storm are less sensitive to the three selected PBL schemes. The main difference, however, is the timing (Figure S20), with simulations using the YSU and QNSE PBL schemes initiating the convective line about 30 min–1 h earlier than the simulation using ACM2. The hybrid nonlocal and local mixing in ACM2 produced a PBL height higher than in YSU and QNSE, which was likely responsible for the later convective initiation. YSU and QNSE PBL schemes initiated convection 1.5 h later than the observed storm, with ACM2 initiating convection about 2.5 h later than observed.

For the chemical mechanisms there are a few differences between the three realizations, due to the inclusion of cloud‐aerosol interactions. Other chemistry interactions (e.g., aerosol‐radiation, cloud chemistry) were turned off. With those options included, there would be increased potential for measurable differences (Figure S17). In Figure , the only distinguishing feature among the three chemical mechanism sensitivity simulations is the shape of the discrete cell to the north of the convective line which is likely due to aerosol‐cloud interactions.

Vertical Extent of Simulated Storms

To assess the sensitivity of storm top altitudes to model parameterization, simulated cloud top and 10 dBZ radar reflectivity echo top box plots are compared to the observed NEXRAD composite 10 dBZ echo tops (Figure ). Cloud tops are determined in WRF simulations as the highest altitude in a column where the cloud‐mixing ratio was at least 0.1 g/kg. Note that the vertical sampling from NEXRAD (1 km) is coarser than the vertical resolution in WRF (250 m). Box plots of WRF output were generated with a coarser resolution to match the NEXRAD resolution and showed no significant difference (not shown). Therefore, it should be noted that WRF results are good to within 1 km of the NEXRAD data.

Overall, the simulated echo tops are underestimated for all BMPs, though the median echo top in MY is slightly closer to that observed. Similar to the BMPs, model‐simulated echo tops are low for all PBL schemes. The echo top distribution with the ACM2 PBL scheme is closer to the observed echo top distribution than the other PBL schemes. Among the chemical mechanisms, there is essentially no difference between them.

Simulated cloud top altitudes are all higher than the observed echo top, which is expected since NEXRAD WSR‐88D radars detect only precipitation‐sized particles. This result is in agreement with previous studies [e.g., Homeyer et al., ; Homeyer, ] of simulated tropopause‐penetrating convection, which show that model simulated echo tops are typically lower than observed, while cloud top altitudes are higher. In contrast to the model simulated echo tops, which showed little variability in the upper bounds (75% and maximum) of the echo top distribution, there are notable differences between the BMPs. NSSL has the highest cloud top altitudes, with maxima more than 1 km above those in MOR.

Chemical Distributions of Simulated Storms

To evaluate the model simulated chemical distribution, simulated vertical profiles of several trace gases superimposed on the DC3 aircraft observations are presented. Vertical profiles were created for the minimum, median, and maximum concentrations over a subsection of the nested domain. The selected area is comparable to the area observed by the aircraft and at comparable times in the storm's evolution (between the observed and simulated storm). In order to analyze only background and anvil‐cloud points, profiles of trace gases were first filtered to remove any convective influence by removing columns where the surface precipitation accumulation is greater than 0 mm and updraft speed is greater than 5 m/s. Next, profiles were stratified into in‐cloud and out‐of‐cloud populations, where in‐cloud are simply those where cloud particles exist in both simulations and observations. The threshold for cloud particles is the same as described in the above section (e.g., cloud indicator for DC8, 2D‐C measurements with threshold of 0.1 g/kg for GV, and total cloud mass mixing ratio with threshold of 0.1 g/kg for WRF simulations).

For water vapor, the median profiles are in general agreement among the microphysics schemes (Figure ), but the maximum mixing ratios differ, especially above the tropopause. Since water vapor injection is an important aspect of STE from extratropical convection due to its chemical and radiative impacts, accurately simulating this process is a key criterion for the evaluation of these simulations. In both in‐cloud and out‐of‐cloud profiles, the MOR and MY schemes show maximum concentrations less than 200 ppmv, while the NSSL scheme shows a maximum concentration as high as 300 ppmv. In comparison, the maximum water vapor concentration observed in the stratosphere by DC3 aircraft was 250 ± 12 ppmv. Thus, it can be concluded that the MOR and MY schemes underrepresent water vapor injection while the NSSL scheme simulates concentrations that were slightly higher than observed.

As noted earlier, the NSSL scheme's design allows for more complex differential size sorting and sedimentation rates compared to the other two BMPs. In particular, simulations with the NSSL BMP had a higher mass mixing ratio of cloud ice in the thunderstorm anvil above the tropopause (Figure ), which likely produced more rapid sublimation and enhancement of water vapor in the stratosphere. Larger frozen hydrometeors (i.e., snow, graupel, and hail) could potentially hydrate the stratosphere, but such hydrometeors are typically located in the interior of the cloud where they require strong upward motion to remain suspended. Ice particles are located near the exterior of the cloud and are directly exposed to the dry stratospheric air, thereby permitting rapid sublimation. MY had the lowest ice mass mixing ratio, but the most snow. MOR BMP has the highest mass mixing ratio of graupel/hail compared to the other two BMPs. These differences in hydrometeor class partitioning are likely responsible for the identified differences in H₂O enhancements in the stratosphere (Figure ).

For the O₃ and CO out‐of‐cloud profiles (Figure ), differences are found between the simulations in the extremes, especially in the depth of mixing above and below the tropopause. The MOR and MY schemes exhibit a similar depth of mixing (3–4 km), while the NSSL scheme shows a shallower layer (1–2 km). It is currently hypothesized that this difference in mixing layer depth is due to the depth of their respective anvil clouds, which extend above the tropopause level. Cross sections, arbitrarily taken, demonstrate the differences in the depths of the simulated storm anvils (Figure ). MOR and MY have similar anvil depths, while NSSL has a shallower anvil. The observations, though somewhat limited, seem to suggest a shallower layer of mixing (~2 km), in agreement with the NSSL scheme.

These simulations also demonstrate the impact of different transport processes (e.g., air mass mixing versus convective injection of ice) as shown in the differences in the peak altitude of the CO/O₃ and H₂O perturbations. Although the aircraft only sampled air ~1 km above the tropopause, WRF‐Chem simulations suggest that H₂O enhancements extended up to 5 km above the tropopause. Transport of constituents into the stratosphere is deeper for H₂O compared to CO and O₃ in part because the cloud and tracer boundaries are often not coincident in WRF‐Chem. More importantly, H₂O enhancements in the stratosphere are sourced by two processes: (1) air mass transport from the troposphere (the only process relevant for O₃ and CO) and (2) rapid sublimation of convectively lofted ice, which can be lofted to higher altitudes than that achieved by convective ascent when gravity wave breaking occurs. In addition, since air is often near saturation and cooling as it ascends, H₂O enhancement from air mass transport is limited to the saturation vapor pressure, which typically reaches 5–10 ppmv at the tropopause. Thus, most of the stratospheric water vapor enhancement is due to the detrainment and sublimation of cloud ice rather than air mass transport. This means that changes in CO and O₃ in the UTLS are largely controlled by a different process than those affecting H₂O, which enables their influence on UTLS composition to vary considerably. More high‐resolution studies are needed to adequately determine this.

For other trace gases (e.g., soluble species such as HNO₃, SO₂, and CH₂O), there is more variability among the BMPs (Figure ). For HNO₃, all three BMPs underestimate the observed maximum HNO₃ concentration in the UTLS by about 1000 pptv for both in‐cloud and out‐of‐cloud profiles, suggesting that the model is over‐scavenging HNO₃, a result consistent with Bela et al. [] and/or underpredicting the lightning‐NO_x production, which subsequently makes HNO₃ via gas‐phase chemistry. MOR and MY have higher HNO₃ concentrations below the tropopause compared to NSSL (and observations), possibly due to the higher precipitation mixing ratios in MOR and MY. The wet scavenging schemes deplete trace gas concentrations based on production of precipitation and therefore more HNO₃ should be depleted at low levels to around 2 km below the tropopause. However, due to the complexities of the cloud physics, it is difficult to say why HNO₃ concentrations are higher in MOR and MY. Explaining these differences should be a part of future studies.

For SO₂, NSSL underestimates concentrations in the UTLS out‐of‐cloud while MOR and MY do a good job reproducing the SO₂ out‐of cloud profile. However, MOR greatly overestimates SO₂ concentrations within cloud. For CH₂O, all three BMPs are in good agreement with observations within cloud. For out‐of‐cloud profiles, MOR and MY overestimate CH₂O while NSSL underestimates CH₂O in the UTLS. Comparison of soluble trace gas profiles within a region not affected by the storm (not shown) reveals very few differences between the BMPs, suggesting that most of the differences (i.e., sensitivity) outlined here are due to interactions with the wet scavenging method employed (see Table ). A deeper investigation into the differences between the parameterizations used is beyond the scope of this study.

Figure  shows the chemical profiles of O₃, H₂O, and CO for the three PBL schemes. There is little impact of PBL scheme choice on the vertical distribution of these trace gases. There are some small differences for CO and other trace gases (e.g., SO₂ and CH₂O (Figure )), but in comparison with the aircraft observations, there is no objectively superior PBL scheme.

Vertical profiles of trace gases exhibit some sensitivity to the choice of chemical mechanism. Overall, there is little sensitivity for O₃, H₂O, and CO (Figure ); however, for other trace gases (Figure ), MOZCART and CBMZ remain consistent, while RACM‐ESRL is measurably different.

In general, MOZCART and CBMZ best reproduce the chemical profiles of HNO₃ and CH₂O while greatly overestimating the maximum SO₂ concentrations in the UTLS. RACM‐ESRL tends to underestimate concentrations of HNO₃ and SO₂ throughout the model domain. In RACM‐ESRL, the lower SO₂ concentrations are due to a SO₂ conversion to SO₄ from cloud chemistry that is not done in MOZCART or CBMZ. The differences in HNO₃, however, are not clear at this time, but potentially due to differences in wet scavenging. An in‐depth look into these differences is the subject of future work.

Discussion

Overall, we have found that WRF‐Chem simulations of the physical characteristics and STE of tropopause‐penetrating convection are most sensitive to the choice of bulk microphysics parameterization. This is not surprising since microphysics play a large role in the structure, organization, and intensity of storms and their ability to hydrate or otherwise modify the composition of the UTLS. Based on the available aircraft observations from the case presented as well as the other two DC3 cases simulated (see supporting information Figures S1–S16), the NSSL 2‐moment BMP best represented the physical and chemical characteristics of the simulated storms. This conclusion was based largely on the scheme's representation of stratospheric water vapor enhancements (from air mass transport and ice injection), which is one of the most important processes that simulations of tropopause‐penetrating convection seek to reproduce. Though not shown, the cloud boundary and O₃/CO boundaries in overshooting convection are often not coincident in WRF‐Chem, which leads to some of the differences in the depth of change in their stratospheric concentrations compared to that for H₂O. Further work should determine whether or not this offset between the cloud and tracer boundaries is representative of reality.

The choice of PBL parameterization shows less sensitivity than the BMP but does show some differences in convective initiation and chemical composition due to differences in mixing and the height of the PBL—which may impact the chemical composition of the air ingested by the convection near cloud base. Among the three PBL schemes, while the physical characteristics were fairly similar, the timing of convective initiation with ACM2 occurred nearly 1 h later than the other two simulations and 2.5 h later than that observed. For the chemical distributions, there were more noticeable, albeit small, differences for the in‐cloud profiles than out‐of‐cloud profiles, though all PBL schemes were very similar. The differences were likely due to differing representations of vertical mixing, which can lead to differences in PBL height and, consequently, air source. Vertical mixing in QNSE is achieved by a local mixing scheme, which gives a lower PBL height, consistent with previous studies comparing local and nonlocal schemes [e.g., Cohen et al., ; Coniglio et al., ]. This could have resulted in a source of air that was less polluted (e.g., lower CO) than the other two schemes (not shown).

While the simulations presented in this study show that the choice in chemical mechanism offers little sensitivity, it should be reiterated that the goals of this study were to assess the sensitivity to bulk transport of mainly passive tracers (e.g., O₃ and CO) and H₂O. Thus, for alternative research questions (e.g., wet scavenging and cloud chemistry) it should be noted that the sensitivity may be larger to the choice of chemical mechanism. Other complexities of the chemical mechanisms, such as gas‐aerosol schemes, would arguably be better suited for studies in which an accurate representation of aerosols is needed. For example, RACM‐ESRL is coupled to the MADE/VBS aerosol module that was designed to give more representative SOA and PM_2.5 concentrations than the other two chemical mechanisms presented here. However, in line with the goals of this study, aerosol‐radiation effects and cloud‐chemistry were turned off. With those interactions enabled, the physical and chemical characteristics would likely have differed, suggesting that the inclusion of these effects offers more sensitivity to STE from extratropical convection than the chosen gas‐phase chemical mechanism (Text S2 and Figure S17). As noted, prognostic aerosols were turned on, though only the RACM‐ESRL and CBMZ chemical mechanism allowed aerosol‐cloud interactions. Turning prognostic aerosols off resulted in higher water vapor perturbations in the stratosphere and a shallower layer of mixing of CO/O₃ out of cloud, and a slightly deeper layer of mixing in‐cloud (Text S3 and Figure S18). For soluble trace gases, there was little difference outside the area of convective influence. Within cloud, the inclusion of aerosol‐cloud interactions resulted in an increase in SO₂ in the UTLS, however, for HNO₃ and HCHO, there were increases in the concentrations of those gases, but the magnitude may vary by chemical mechanism (Text S3 and Figure S19).

Conclusions

Of the three different model parameterization types tested in this study, WRF‐Chem simulations of tropopause‐penetrating convection and STE are most sensitive to the choice in BMP. As expected, there is measurable sensitivity of the organization and vertical extent of simulated convection to the choice of BMP. Furthermore, the simulations in this study show that convectively injected water into the stratosphere is also sensitive to the choice of BMP. Among the three BMPs tested, the NSSL 2‐moment scheme provided the best representation of both the observed physical characteristics of the storm and the composition of the UTLS.

There was little sensitivity of the physical characteristics of the storms to the chosen chemical mechanism. For assessing bulk transport of mostly passive trace gases, there were little apparent differences. Thus, the choice in chemical mechanism should be based on the specific research questions the user wants to investigate. For the PBL schemes tested, there was little sensitivity in both the physical structure of the simulated storm and composition of the UTLS, though there were small differences in the timing of convective initiation.

Acknowledgments

We thank Ted Mansell and the two anonymous reviewers for their helpful suggestions in improving the original manuscript. The authors acknowledge ECMWF for providing the ERA‐Interim reanalysis output, which were obtained from the Research Data Archive (RDA) maintained by the Computational and Informational Systems Laboratory (CISL) at the National Center for Atmospheric Research (NCAR); the original data are available at http://has.ncdc.noaa.gov/ and http://rda.ucar.edu/, respectively. We also acknowledge use of MOZART‐4 global model output available at http://www.acom.ucar.edu/wrf‐chem/mozart.shtml, the use of the WRF‐Chem preprocessor tool {mozbc, fire_emiss, etc.} provided by the Atmospheric Chemistry Observations and Modeling Lab (ACOM) of NCAR, and Stu Mckeen for (NOAA/ESRL) for providing 4 km gridded anthropogenic emissions from the EPA Inventory, and the National Center for Atmospheric Research is supported by the National Science Foundation. The computing for this project was performed at the OU Supercomputing Center for Education & Research (OSCER) at the University of Oklahoma (OU). This work was supported by the National Science Foundation grant AGS‐1522910.