7. Extensions to CAM¶

7.1. Introduction¶

This section contains a description of the neutral constituent chemical options available CAM and WACCM4.0, including different chemical schemes, emissions, boundary conditions, lightning, dry depositions and wet removal; 2) the photolysis approach; 3) numerical algorithms used to solve the corresponding set of ordinary differential equations.; 4) additions to superfast chemistry.

7.2. Chemistry¶

7.2.1. Chemistry Schemes¶

For CAM-Chem, an extensive tropospheric chemistry option is available (trop mozart), as well as an extensive tropospheric and stratospheric chemistry (trop-strat mozart) as discussed in detail in ([LEH+12]), including a list of all species and reactions. Furthermore, a superfast chemistry option is available for CAM, as discussed in Section Superfast Chemistry. For each chemical scheme, CAM-chem uses the same chemical preprocessor as MOZART-4. This preprocessor generates Fortran code for each specific chemical mechanism, allowing for an easy update and modification of existing chemical mechanisms. In particular, the generated code provides two chemical solvers, one explicit and one semi-implicit, which the user specifies based on the chemical lifetime of each species. For all supported compsets, this generated code is available in a sub-directory of atm/src/chemistry.

7.2.1.1. The Bulk Aerosol Model¶

CAM4-chem uses the bulk aerosol model discussed in [LHE+05] and [EWH+10]. This model has a representation of aerosols based on the work by [{tie:02] and [TMW+05a], i.e. sulfate aerosol is formed by the oxidation of SO $_{2}$ in the gas phase (by reaction with the hydroxyl radical) and in the aqueous phase (by reaction with ozone and hydrogen peroxide). Furthermore, the model includes a representation of ammonium nitrate that is dependent on the amount of sulfate present in the air mass following the parameterization of gas/aerosol partitioning by [MDPL02]. Because only the bulk mass is calculated, a lognormal distribution is assumed for all aerosols using different mean radius and geometric standard deviation [LAC+03]. The conversion of carbonaceous aerosols (organic and black) from hydrophobic to hydrophilic is assumed to occur within a fixed 1.6 days [TMW+05a]. Natural aerosols (desert dust and sea salt) are implemented following [MLTW06] and [MML+06] and the sources of these aerosols are derived based on the model calculated wind speed and surface conditions. In addition, secondary-organic aerosols (SOA) are linked to the gas-phase chemistry through the oxidation of atmospheric non-methane hydrocarbons (NMHCs), as in [LTB+04].

7.2.1.3. Trop MOZART Chemistry¶

The extensive tropospheric chemistry scheme represents a minor update to the MOZART-4 mechanism, fully described in ([EWH+10]). In particular, we have included C $_{2}$ H $_{2}$ , HCOOH, HCN and CH $_{3}$ CN. Reaction rates have been updated to JPL-2006 ([SanderSPeal06]). A minor update has been made to the isoprene oxidation scheme, including an increase in the production of glyoxal. This mechanism is mainly of relevance in the troposphere and is intended for simulations for which long-term trends in the stratospheric composition are not crucial. Therefore, in this configuration, the stratospheric distributions of long-lived species (see discussion below) are specified from previously performed WACCM simulations ([GMK+07]); see Section Boundary conditions).

7.2.1.4. Trop-Strat MOZART Chemistry¶

The extensive tropospheric and stratospheric chemistry includes the full stratospheric chemistry from WACCM4.0, with an updated enforcement of the conservation of total chlorine and total bromine under advection to improve the performance of the model in simulating the ozone hole. In addition, we have updated the heterogeneous chemistry module to reflect that the model was underestimating the supercooled ternary solution (STS) surface area density (SAD), see more detail in Section [sec:waccm]; (???), Kinnison et al, 2012, (in prepareation).

7.2.1.5. SOA calculation in CAM-Chem¶

An SOA simulation of intermediate complexity is also available in CAM-Chem. This is based on the 2-product model scheme of [OJG+97], as implemented in CAM-Chem by [HHH+08]. This treats the products of VOC oxidation as semi-volatile species, which re-partition every time step based on the temperature (enthalpy of vaporization of 42 kJmol-1) and organic aerosol mass available for condensation of vapours [Pan94]. In CAM-Chem we treat secondary organic aerosol formation from the products of isoprene, monoterpene and aromatic (benzene, toluene and xylene) oxidation by OH, O $_{3}$ and NO $_{3}$ . The yields and partitioning coefficients are based on smog chamber studies [GCFS99][HSN+08][NKC+07]. The SOA calculation is setup to run with biogenic emissions calucated by the MEGAN2.1 model (see Section Emissions).

7.2.2. Emissions¶

Surface emissions are used in as a flux boundary condition for the diffusion equation of all applicable tracers in the planetary boundary-layer scheme. The surface flux files used in the released version are discussed in (???) and conservatively remapped from their original resolution (monthly data available every decade at 0.5x0.5) to (monthly data every year at 1.9x2.5). The remapping is made offline to avoid the internal remapping, which consists of a simple linear interpolation and therefore does not ensure exact conservation of emissions between resolutions.

7.2.2.1. Emissions of Trace Gases¶

Emissions for historic and future model simulations are based on ACCMIP ((???)) and different RCP scenarios, which are available for the years 1850-2000, and 2000-2100.

Additional emissions are available for a short period covering 1992-2010, as discussed in (???). More specifically, for 1992-1996, which is prior to satellite-based fire inventories, monthly mean averages of the fire emissions for 1997-2008 from GFED2 (??? and updates) are used for each year. For 2009-2010, fire emissions are from FINN (Fire INventory from NCAR) (???). If running with FINN fire emissions, additional species are availabel: NO $_{2}$ , BIGALD, CH $_{3}$ COCHO, CH $_{3}$ COOH, CRESOL, GLYALD, HYAC, MACR, MVK. Most of the anthropogenic emissions come from the POET (Precursors of Ozone and their Effects in the Troposphere) database for 2000 (???). The anthropogenic emissions (from fossil fuel and biofuel combustion) of black and organic carbon determined for 1996 are from (???). For SO $_{2}$ and NH $_{3}$ , anthropogenic emissions are from the EDGAR-FT2000 and EDGAR-2 databases, respectively (http://www.mnp.nl/edgar/).

For Asia, these inventories have been replaced by the Regional Emission inventory for Asia (REAS) with the corresponding annual inventory for each year simulated (???). Only the Asian emissions from REAS vary each year, all other emissions are repeated annually for each year of simulation. The DMS emissions are monthly means from the marine biogeochemistry model HAMOCC5, representative of the year 2000 (???).

Additional emissions (volcanoes and aircraft) are included as three-dimensional arrays, conservatively-remapped to the CAM-chem grid. The volcanic emission are from (???) and the aircraft (NO:math:_{2}) emissions are from (???). In the case of volcanic emissions (SO:math:_{2} and SO $_{2}$ ), an assumed 2% of the total sulfur mass is directly released as SO $_{2}$ . SO $_{2}$ emissions from continuously outgassing volcanoes are from the GEIAv1 inventory (Andres and Kasgnoc, 1998). Totals for each year and emitted species are listed in (???), Table 7. Aerosol Emissions available to be used in CAM5-Chem are described above (Section 4.8.1.).

7.2.2.2. Biogenic emissions¶

Biogenic emissions can be calculated by the Model of Emissions of Gases and Aerosols from Nature version 2.1 (MEGAN2.1) (???). In this case, MEGAN2.1 is coupled to the CESM atmosphere and land model. Biogenic emissions of volatile organic compounds (i.e. isoprene and monoterpenes) are calculated based upon emission factors, land cover (LAI and PFT), and driving meteorological variables. CO $_{2}$ effect on isoprene emission is also included (???). Emission factors of different MEGAN compounds can be specified from mapped files or based on PFTs. These are made available for atmospheric chemistry, unless the user decides to explicitly set those emissions using pre-defined (i.e. contained in a file) gridded values. Details of this implementation in the CLM3 are discussed in (???) and (???): Vegetation in the CLM model is described by 17 plant function types (PFTs, see (??? Table 1)). Present-day land cover data such as leaf area index are consistent with MODIS land surface data sets (???). Alternate land cover and density can be either specified or interactively simulated with the dynamic vegetation model (CLMCNDV) or the carbon nitrogen model (CLMCN) of the CLM for any time period of interest. Additional namelist parameters have been included to facilitate the mapping between the emissions in MEGAN2.1 (147 species) and the chemical mechanism. Surface emissions without biogenic emissions have to be used if the MEGAN2.1 model produces biogenic emissions to prevent double counting.

7.2.3. Boundary conditions¶

7.2.3.1. Lower boundary conditions¶

For all long-lived species (methane and longer lifetimes, in addition to hydrogen and methyl bromide) (??? see Table 3), the surface concentrations are specified using the historical reconstruction from (???). In addition, for CO $_{2}$ and CH $_{4}$ , an observationally-based seasonal cycle and latitudinal gradient are imposed on the annual average values provided by (???). These values are used in the model by overwriting at each time step the corresponding model mixing ratio in the lowest model level with the time (and latitude, if applicable) interpolated specified mixing ratio.

7.2.3.2. Specified stratospheric distributions¶

For the trop-mozart chemistry, no stratospheric chemistry is explicitly represented in the model. Therefore, it is necessary to ensure a proper distribution of some chemically-active stratospheric (namely O $_{3}$ , NO, NO $_{2}$ , HNO $_{3}$ , CO, CH $_{4}$ , N $_{2}$ O, and N $_{2}$ O $_{5}$ ) species, as is the case for MOZART-4. This monthly-mean climatological distribution is obtained from WACCM simulations covering 1950-2005 (???). Because of the vast changes that occur over that time period, our data distribution provides files for three separate periods: 1950-1959, 1980-1989 and 1996-2005. This ensure that users can perform simulations with a stratospheric climatology representative of the pre-CFC era, as well as during the high CFC and post-Pinatubo era. Additional datasets for different RCP runs are also available or can easily be constructed if necessary.

7.2.3.3. Upper boundary condition¶

The model top at about 40km is considered a rigid lid (no flux across that boundary) for all chemical species. For trop-mozart

7.2.4. Lightning¶

The emissions of NO from lightning are included as in (???), i.e. using the Price parameterization ((???; ???), scaled to provide a global annual emission of 3-4 Tg(N)/year. The vertical distribution follows (???) as in (???). In addition, the strength of intra-cloud (IC) lightning strikes is assumed to be equal to cloud-to-ground strikes, as recommended by (???).

Lightning NOx can be modifid in the namelist. For CAM5-Chem, lightning NOx is increased by a factor of 3 to reach the same emissions of 3-4 Tg(N)/year.

7.2.5. Dry deposition¶

Dry deposition is represented following the resistance approach originally described in (???), as discussed in (???), this earlier paper was subsequently updated and we have included all updates (???; ???). Following this approach, all deposited chemical species (the specific list of deposited species is defined along with the chemical mechanisms, see section 4) are mapped to a weighted-combination of ozone and sulfur dioxide depositions; this combination represents a definition of the ability of each considered species to oxidize or to be taken up by water. In particular, the latter is dependent on the effective Henry’s law coefficient. While this weighting is applicable to many species, we have included specific representations for CO/H $_{2}$ (???; ???) and peroxyacetylnitrate (PAN) (???). Furthermore, it is assumed that the surface resistance for SO $_{2}$ can be neglected (???). Finally, following (???), the deposition velocities of black and organic carbonaceous aerosols are specified to be 0.1 cm/s over all surfaces. Dust and sea-salt are represented following (???) and (???).

The computation of deposition velocities in CAM-chem takes advantage of its coupling to the Community Land Model (CLM; http://www.cesm.ucar.edu/ models/cesm1.0/ clm/index.shtml). In particular, the computation of surface resistances in CLM leads to a representation at the level of each plant functional type (Table 1) of the various drivers for deposition velocities. The grid-averaged velocity is computed as the weighted-mean over all land cover types available at each grid point. This ensures that the impact on deposition velocities from changes in land cover, land use or climate is taken into account. All species in the mechanism are per default affected by dry deposition if depostion velocities are defined in the model.

7.2.6. Wet removal¶

Wet removal of soluble gas-phase species is the combination of two processes: in-cloud, or nucleation scavenging (rainout), which is the local uptake of soluble gases and aerosols by the formation of initial cloud droplets and their conversion to precipitation, and below-cloud, or impaction scavenging (washout), which is the collection of soluble species from the interstitial air by falling droplets or from the liquid phase via accretion processes (e.g. ???).

Removal is modeled as a simple first-order loss process X $_{iscav}$ =X $_{i} \cdot$ F $\cdot (1- exp(- \lambda ~ \Delta$ t)). In this formula, X $_{iscav}$ is the species mass (in kg) and X $_{i}$ scavenged in time step $\Delta$ t , F is the fraction of the grid box from which tracer is being removed, and $\lambda$ is the loss rate. In-cloud scavenging is proportional to the amount of condensate converted to precipitation, and the loss rate depends on the amount of cloud water, the rate of precipitation formation, and the rate of tracer uptake by the liquid phase. Below-cloud scavenging is proportional to the precipitation flux in each layer and the loss rate depends on the precipitation rate and either the rate of tracer uptake by the liquid phase (for accretion processes), the mass-transfer rate (for highly soluble gases and small aerosols), or the collision rate (for larger aerosols).

In CAM-chem two separate parameterizations are available: (???) from MOZART-2 and (???). The distinguishing features of the Neu and Prather scheme are related to three aspects of the parameterization: 1) the partitioning between in-cloud and below cloud scavenging, 2) the treatment of soluble gas uptake by ice and 3) the Neu and Prather scheme uniquely accounts for the spatial distribution of clouds in a column and the overlap of condensate and precipitation. Given a cloud fraction and precipitation rate in each layer, the scheme determines the fraction of the gridbox exposed to precipitation from above and that exposed to new precipitation formation under the assumption of maximum overlap of the precipitating fraction. Each model level is partitioned into as many as four sections, each with a gridbox fraction, precipitation rate, and precipitation diameter: 1) Cloudy with precipitation falling through from above; 2) Cloudy with no precipitation falling through from above; 3) Clear sky with precipitation falling through from above; 4) Clear sky with no precipitation falling from above. Any new precipitation formation is spread evenly between the cloudy fractions (1 and 2). In region 3, we assume a constant rate of evaporation that reduces both the precipitation area and amount so that the rain rate remains constant. Between levels, we average the properties of the precipitation and retain only two categories, precipitation falling into cloud and precipitation falling into ambient air, at the top boundary of each level. If the precipitation rate drops to zero, we assume full evaporation and random overlap with any precipitating levels below. Our partitioning of each level and overlap assumptions are in many ways similar to those used for the moist physics in the ECMWF model (???).

The transfer of soluble gases into liquid condensate is calculated using Henry’s Law, assuming equilibrium between the gas and liquid phase. Nucleation scavenging by ice, however, is treated as a burial process in which trace gas species deposit on the surface along with water vapor and are buried as the ice crystal grows. (???) have found that the burial model successfully reproduces the molar ratio of HNO $_{3}$ to H $_{2}$ O on ice crystals as a function of temperature for a large number of aircraft campaigns spanning a wide variety of meteorological conditions. We use the empirical relationship between the HNO $_{3}$ H $_{2}$ O molar ratio and temperature given by (???) to determine in-cloud scavenging during ice particle formation, which is applied to nitric acid only. Below-cloud scavenging by ice is calculated using a rough representation of the riming process modeled as a collision-limited first order loss process. (???) provide a full description of the scavenging algorithm.

On the other hand, the Horowitz approach uses the rain generation diagnostics from the large-scale and convection precipitation parameterizations in CAM; equilibrium between gas-phase and liquid phase is then assumed based on the effective Henry’s law.

7.3. Photolytic Approach (Neutral Species)¶

The calculation of the photolysis coefficients is divided into two regions: (1) 120 nm to 200 nm (33 wavelength intervals); (2) 200 nm to 750 nm (67 wavelength intervals). The total photolytic rate constant (J) for each absorbing species is derived during model execution by integrating the product of the wavelength dependent exo-atmospheric flux (F:math:_{exo}); the atmospheric transmission function (or normalized actinic flux) (N:math:_A), which is unity at the top of atmosphere in most wavelength regions; the molecular absorption cross-section ( $\sigma$ ); and the quantum yield ( $\phi$ ). The exo-atmospheric flux over these wavelength intervals can be specified from observations and varied over the 11-year solar sunspot cycle (see section [sec:short_wave]).

The wavelength-dependent transmission function is derived as a function of the model abundance of ozone and molecular oxygen. For wavelengths greater than 200 nm a normalized flux lookup table (LUT) approach is used, based on the 4-stream version of the Stratosphere, Troposphere, Ultraviolet (STUV) radiative transfer model (S. Madronich, personal communication), (???). The transmission function is interpolated from the LUT as a function of altitude, column ozone, surface albedo, and zenith angle. The temperature and pressure dependences of the molecular cross sections and quantum yields for each photolytic process are also represented by a LUT in this wavelength region. At wavelengths less than 200 nm, the wavelength-dependent cross section and quantum yields for each species are specified and the transmission function is calculated explicitly for each wavelength interval. There are two exceptions to this approach. In the case of J(NO) and J(O $_2$ ), detailed photolysis parameterizations are included inline. In the Schumann-Runge Band region (SRBs), the parameterization of NO photolysis in the $\delta$ -bands is based on (???). This parameterization includes the effect of self-absorption and subsequent attenuation of atmospheric transmission by the model-derived NO concentration. For J(O $_2$ ), the SRB and Lyman-alpha parameterizations are based on (???) and (???), respectively.

While the lookup table provides explicit quantum yields and cross-sections for a large number of photolysis rate determinations, additional ones are available by scaling of any of the explicitly defined rates. This process is available in the definition of the chemical preprocessor input files (see (??? Table 3) for a complete list of the photolysis rates available). The impact of clouds on photolysis rates is parameterized following (???). However, because we use a lookup table approach, the impact of aerosols (tropospheric or stratospheric) on photolysis rates cannot be represented.

As an extension of MOZART-4 and to provide the ability to seamlessly perform tropospheric and stratospheric chemistry simulations, the calculation of photolysis rates for wavelengths shorter than 200 nm is included; this was shown to be important for ozone chemistry in the tropical upper troposphere (???). In addition, because the standard configuration of CAM only extends into the lower stratosphere (model top is usually around 40 km), an additional layer of ozone and oxygen above the model top is included to provide a very accurate representation of photolysis rates in the upper portion of the model as compared to the equivalent calculation using a fully-resolved stratospheric distribution.

In addition, tropospheric photolysis rates can be computed interactively. Users interested in using this capability have to contact the Chemistry-CLimate Working Group Liaison as this is an unsupported option.

7.4. Numerical Solution Approach¶

Chemical and photochemical processes are expressed by a system of time-dependent ordinary differential equations at each point in the spatial grid, of the following form:

$\frac{d\vec{y}}{dt} = \vec{P}(\vec{y}, t) - \vec{L}(\vec{y}, t) \cdot \vec{y}$

$\vec y(t) = \{y_i(t)\} \quad i = 1, 2, \ldots, N$

where $\vec y$ is the vector of all solution variables (chemical species), $N$ is the number of variables in the system, and $y_i$ represents the $i^{th}$ variable. $\vec P$ and $\vec L$ represent the production and loss rates, which are, in general, non-linear functions of the $y_i$ . This system of equations is solved via two algorithms: an explicit forward Euler method:

$y_i^{n+1} = y_i^n + \Delta t \cdot f_i(t_{n}, y^{n})$

in the case of species with long lifetimes and weak forcing terms (e.g., N $_2$ O), and a more robust implicit backward Euler method:

$y_i^{n+1} = y_i^n + \Delta t\cdot f_i(t_{n+1}, y^{n+1})$

for species that comprise a``stiff system” with short lifetimes and strong forcings (e.g., OH). Here $n$ represents the time step index. Each method is first order accurate in time and conservative. The overall chemistry time step, $\Delta t = t_{n+1}-t_n$ , is fixed at 30 minutes. Preprocessing software requires the user to assign each solution variable, $y_i$ , to one of the solution schemes. The discrete analogue for methods (2) and (3) above results in two systems of algebraic equations at each grid point. The solution to these algebraic systems for equation (2) is straightforward (i.e., explicit). The algebraic system from the implicit method (3) is quadratically non-linear. This system can be written as:

$\vec{G}(\vec{y}^{\,\, n+1})=\vec{y}^{\,\, n+1}-\vec{y}^{\,\, n}- \Delta t\cdot\vec{f}(t_{n+1},\vec{y}^{\,\, n+1})=0$

Here $G$ is an $N$ -valued, non-linear vector function, where $N$ equals the number of species solved via the implicit method. The solution to equation (4) is solved with a Newton- Raphson iteration approach as shown below:

$\vec{y}^{\,\, n+1}_{m+1} = \vec{y}^{\,\, n+1}_m - \vec{J} \cdot \vec{G}(\vec{y}^{\,\, n+1}_m); \; m=0,1,\ldots, M$

Where $m$ is the iteration index and has a maximum value of ten. The elements of the Jacobian matrix $\vec J$ are given by:

$J_{ij}=\frac{\partial G_i}{\partial y_j}$

The iteration and solution of equation (5) is carried out with a sparse matrix solution algorithm. This process is terminated when the given solution variable changes in a relative measure by less than a prescribed fractional amount. This relative error criterion is set on a species by species basis, and is typically 0.001; however, for some species (e.g., O $_3$ ), where a tighter error criterion is desired, it is set to 0.0001. If the iteration maximum is reached (for any species) before the error criterion is met, the time step is cut in half and the solution to equation (5) is iterated again. The time step can be reduced five times before the solution is accepted. This approach is based on the work of [SPDIC96] and [SVvL+97]; see also [BOT99].

7.5. Superfast Chemistry¶

7.5.1. Chemical mechanism¶

The super-fast mechanism was developed for long coupled chemistry-climate simulations, and is based on an updated version of the full non-methane hydrocarbon effects (NMHC) chemical mechanism for the troposphere and stratosphere used in the Lawrence Livermore National Laboratory off-line 3D global chemistry-transport model (IMPACT) citep[]rotman:04. The super-fast mechanism includes 15 photochemically active trace species (O:math:_{3}, OH, HO $_{2}$ , H $_{2}$ O $_{2}$ , NO, NO $_{2}$ , HNO $_{3}$ , CO, CH $_{2}$ O, CH $_{3}$ O $_{2}$ , CH $_{3}$ OOH, DMS, SO $_{2}$ , SO $_{4}$ , and C $_{5}$ H $_{8}$ ) that allow us to calculate the major terms by which global change operates in tropospheric ozone and sulfate photochemistry. The families selected are Ox, HOx, NOy, the CH $_{4}$ oxidation suite plus isoprene (to capture the main NMHC effects), and a group of sulfur species to simulate natural and anthropogenic sources leading to sulfate aerosol. Sulfate aerosols is handled following [TMW+05b]. In this scheme, CH4 concentrations are read in from a file and uses CAM3.5 simulations [LamarqueBondEyring+10]. The super-fast mechanism was validated by comparing the super-fast and full mechanisms in side-by-side simulations.

7.5.2. Emissions for CAM4 superfast chemistry¶

|lccc| & Anthro. &Natural & Interactive
CHO & x & x &
CO & x & x &
DMS & & x &
ISOP & & & x
NO & x & &
SO & x & &

7.5.3. LINOZ¶

Linoz is linearized ozone chemistry for stratospheric modeling [MOH+00]. It calculates the net production of ozone (i.e., production minus loss) as a function of only three independent variables: local ozone concentration, temperature, and overhead column ozone). A zonal mean climatology for these three variables as well as the other key chemical variables such a total odd-nitrogen methane abundance is developed from satellite and other in situ observations. A relatively complete photochemical box model [Pra92] is used to integrate the radicals to a steady state balance and then compute the net production of ozone. Small perturbations about the chemical climatology are used to calculate the coefficients of the first-order Taylor series expansion of the net production in terms of local ozone mixing ratio (f), temperature (T), and overhead column ozone (c).

$\begin{aligned} \frac{df}{df} &=& (P - L)^o + \left.{\frac{\delta (P - L)}{\delta f}}\right|_o(f - f^o) + \left.\frac{\delta (P - L)}{\delta T}\right|_o (T - T^o)\\ \nonumber & & + \left.\frac{\delta (P - L)}{\delta c}\right|_o(c - c^o) \end{aligned}$

The photochemical tendency for the climatology is denoted by $(P-L)_o$ , and the climatology values for the independent variables are denoted by $f_o$ , $c_o$ , and $T_o$ , respectively. Including these four climatology values and the three partial derivatives, Linoz is defined by seven tables. Each table is specified by 216 atmospheric profiles: 12 months by 18 latitudes (85:math:^oS to 85 $^o$ N). For each profile, quantities are evaluated at every 2 km in pressure altitude from $z^*$ = 10 to 58 km ( $z^*$ = 16 km log $_10$ (1000/p)). These tables (calculated for each decade, 1850-2000 to take into account changes in CH4 and N2O) are automatically remapped onto the CAM-chem grid with the mean vertical properties for each CAM-chem level calculated as the mass-weighted average of the interpolated Linoz profiles. Equation (1) is implemented for the chemical tendency of ozone in CAM-chem.

7.5.4. Parameterized PSC ozone loss¶

In the superfast chemistry, we incorporate the PSCs parameterization scheme of [CLBRG90] when the temperature falls below 195 K and the sun is above the horizon at stratospheric altitudes. The O $_3$ loss scales as the squared stratospheric chlorine loading (normalized by the 1980 level threshold). In this formulation PSC activation invokes a rapid e-fold of O $_3$ based on a photochemical model, but only when t he temperature stays below the PSC threshold. The stratospheric chlorine loading (1850-2005) is input in the model using equivalent effective stratospheric chlorine (EESC) ([NDWN07]) table based on observed mixing ratios at the surface.

This can be used instead of the more explicit representation available from WACCM in the strat-trop configuration.

7.6. Physical Parameterizations¶

In WACCM4.0, we extend the physical parameterizations used in CAM4 by adding constituent separation velocities to the molecular (vertical) diffusion and modifying the gravity spectrum parameterization. Both of these parameterizations are present, but not used, in CAM4. In addition, we replace the CAM4 parameterizations for both solar and longwave radiation above $\sim 65$ km, and add neutral and ion chemistry models.

7.6.1. Domain and Resolution¶

WACCM4.0 has 66 vertical levels from the ground to $5.1 \times 10^{-6}$ hPa, as in the previous WACCM versions. As in CAM4, the vertical coordinate is purely isobaric above 100 hPa, but is terrain following below that level. At any model grid point, the local pressure p is determined by

$p(i,j,k) = A(k) \, p_0 + B(k) \, p_s(i,j)$

where $A$ and $B$ are functions of model level, $k$ , only; $p_0=10^3$ hPa is a reference surface pressure; and $p_s$ is the predicted surface pressure, which is a function of model longitude and latitude (indexed by $i$ and $j$ ). The finite volume dynamical core uses locally material surfaces for its internal vertical coordinate and remaps (conservatively interpolates) to the hybrid surfaces after each time step.

Within the physical and chemical parameterizations, a local pressure coordinate is used, as described by (6) . However, in the remainder of this note we refer to the vertical coordinate in terms of log-pressure altitude

$Z = H \log\left(\frac{p_0}{p}\right).$

The value adopted for the scale height, $H=7$ km, is representative of the real atmosphere up to $\sim 100$ km, above that altitude temperature increases very rapidly and the typical scale height becomes correspondingly larger. It is important to distinguish $Z$ from the geopotential height $z$ , which is obtained from integration of the hydrostatic equation.

In terms of log-pressure altitude, the model top level is found at $Z=140$ km ( $z\simeq 150$ km). It should be noted that the solution in the top 15-20 km of the model is undoubtedly affected by the presence of the top boundary. However, it should not be thought of as a sponge layer, since molecular diffusion is a real process and is the primary damping on upward propagating waves near the model top. Indeed, this was a major consideration in moving the model top well above the turbopause. Considerable effort has been expended in formulating the upper boundary conditions to obtain realistic solutions near the model top and all of the important physical and chemical processes for that region have been included.

The standard vertical resolution is variable; it is 3.5 km above about 65 km, 1.75 km around the stratopause (50 km), 1.1-1.4 km in the lower stratosphere (below 30 km), and 1.1 km in the troposphere (except near the ground where much higher vertical resolution is used in the planetary boundary layer).

Two standard horizontal resolutions are supported in WACCM4.0: the $4 \times 5 ^{\circ}$ (latitude $\times$ longitude) low resolution version has 72 longitude and 46 latitude points; the $1.9 \times 2.5 ^{\circ}$ medium resolution version has 96 longitude and 144 latitude points. A $0.9 \times 2.5 ^{\circ}$ high resolution version of has had limited testing, and is not yet supported, due to computational cost constraints. The $4 \times 5 ^{\circ}$ version has been used extensively for MLT studies, where it gives very similar results to the $1.9 \times 2.5 ^{\circ}$ version. However, caution should be exercised in using $4 \times 5 ^{\circ}$ results below the stratopause, since the meridional resolution may not be sufficient to represent adequately the dynamics of either the polar vortex or synoptic and planetary waves.

At all resolutions, the time step is 1800 s for the physical parameterizations. Within the finite volume dynamical core, this time step is subdivided as necessary for computational stability.

7.6.2. Molecular Diffusion and Constituent Separation¶

The vertical diffusion parameterization in CAM4 provides the interface to the turbulence parameterization, computes the molecular diffusivities (if necessary) and finally computes the tendencies of the input variables. The diffusion equations are actually solved implicitly, so the tendencies are computed from the difference between the final and initial profiles. In WACCM4.0, we extend this parameterization to include the terms required for the gravitational separation of constituents of differing molecular weights. The formulation for molecular diffusion follows (???)

A general vertical diffusion parameterization can be written in terms of the divergence of diffusive fluxes:

$\begin{aligned} \frac{\partial }{\partial t} (u,v,q) &=& - \frac{1}{\rho} \frac{\partial}{\partial z} (F_u, F_v, F_q) \\ \end{aligned}$

$\begin{aligned} \frac{\partial}{\partial t} s &=& - \frac{1}{\rho} \frac{\partial}{\partial z} F_H + D \end{aligned}$

where $s = c_p T + g z$ is the dry static energy, $z$ is the geopotential height above the local surface (does not include the surface elevation) and $D$ is the heating rate due to the dissipation of resolved kinetic energy in the diffusion process. The diffusive fluxes are defined as:

$F_{u,v} =-\rho K_m \frac{\partial}{\partial z}(u,v),$

$F_{H} =-\rho K_H \frac{\partial s}{\partial z} +\rho K_H^t\gamma_{H} ,$

$F_{q} =-\rho K_q \frac{\partial q}{\partial z} +\rho K_q^t\gamma_{q} + {\rm sep-flux} .$

The viscosity $K_m$ and diffusivities $K_{q,H}$ are the sums of: turbulent components $K_{m,q,H}^t$ , which dominate below the mesopause; and molecular components $K_{m,q,H}^m$ , which dominate above 120 km. The non-local transport terms $\gamma_{q,H}$ are given by the ABL parameterization and and the kinetic energy dissipation is

$D \equiv -\frac{1}{\rho} \left( F_u\frac{\partial u}{\partial z} + F_v\frac{\partial v}{\partial z} \right).$

The treatment of the turbulent diffusivities $K_{m,q,H}^t$ , the energy dissipation $D$ and the nonlocal transport terms $\gamma_{H,q}$ is described in the Technical Description and will be omitted here.

7.6.2.1. Molecular viscosity and diffusivity¶

The empirical formula for the molecular kinematic viscosity is

$K_m^m = 3.55\times 10^{-7} T^{2/3} / \rho,$

and the molecular diffusivity for heat is

$K_H^m = P_r K_m^m,$

where $P_r$ is the Prandtl number and we assume $P_r=1$ in WACCM4.0. The constituent diffusivities are

$K_q^m = T^{1/2} M_w/ \rho,$

where $M_w$ is the molecular weight.

7.6.2.2. Diffusive separation velocities¶

As the mean free path increases, constituents of different molecular weights begin to separate in the vertical. In WACCM4.0, this separation is represented by a separation velocity for each constituent with respect mean air. Since extends only into the lower thermosphere, we avoid the full complexity of the separation problem and represent mean air by the usual dry air mixture used in the lower atmosphere ( $M_w = 28.966$ ) ([BK73]).

7.6.2.3. Discretization of the vertical diffusion equations¶

In CAM4, as in previous version of the CCM, (7) [eq:fluxz2]) are cast in pressure coordinates, using

$dp = -\rho g dz,$

and discretized in a time-split form using an Euler backward time step. Before describing the numerical solution of the diffusion equations, we define a compact notation for the discrete equations. For an arbitrary variable $\psi$ , let a subscript denote a discrete time level, with current step $\psi_n$ and next step $\psi_{n+1}$ . The model has $L$ layers in the vertical, with indexes running from top to bottom. Let $\psi^k$ denote a layer midpoint quantity and let $\psi^{k\pm}$ denote the value at the interface above (below) $k$ . The relevant quantities, used below, are then:

$\begin{aligned} \psi^{k+}&=&(\psi^{k}+\psi^{k+1})/2,\quad k\in (1,2,3,...,L-1) \nonumber \\ \psi^{k-}&=&(\psi^{k-1}+\psi^{k})/2,\quad k\in (2,3,4...,L) \nonumber \\ \delta^{k}\psi &=& \psi^{k+}-\psi^{k-}, \nonumber \\ \delta^{k+}\psi &=& \psi^{k+1}-\psi^{k}, \nonumber \\ \delta^{k-}\psi &=& \psi^{k}-\psi^{k-1}, \nonumber \\ \psi_{n+}&=&(\psi_{n}+\psi_{n+1})/2, \nonumber \\ \delta_n\psi &=& \psi_{n+1}-\psi_{n}, \nonumber \\ \delta t &=& t_{n+1}-t_{n}, \nonumber \\ \Delta^{k,l} &=& 1,\ k=l, \nonumber \\ &=& 0,\ k\neq l. \nonumber\end{aligned}$

Like the continuous equations, the discrete equations are required to conserve momentum, total energy and constituents. Neglecting the nonlocal transport terms, the discrete forms of (7) [eq:contz3]) are:

$\frac{\delta_n (u,v,q)^k}{\delta t} = g \frac{\delta^k F_{u,v,q}}{\delta^k p}$

$\frac{\delta_n s^k}{\delta t} = g \frac{\delta^k F_{H}}{\delta^k p} + D^k.$

For interior interfaces, $1\le k \le L-1$ ,

$F_{u,v}^{k+} = \left(g\rho^2 K_m\right)_n^{k+} \frac{\delta^{k+} (u,v)_{n+1}} {\delta^{k+} p}$

$F_{q,H}^{k+} = \left(g\rho^2 K_{q,H}\right)_n^{k+} \frac{\delta^{k+} (u,v)_{n+1}} {\delta^{k+} p}.$

Surface fluxes $F_{u,v,q,H}^{L+}$ are provided explicitly at time $n$ by separate surface models for land, ocean, and sea ice while the top boundary fluxes are usually $F_{u,v,q,H}^{1-}=0$ . The turbulent diffusion coefficients $K_{m,q,H}^{t}$ and non-local transport terms $\gamma_{q,H}$ are calculated for time $n$ by the turbulence model (identical to CAM4). The molecular diffusion coefficients, given by (13) [Kmq]) are also evaluated at time $n$ .

7.6.2.4. Solution of the vertical diffusion equations¶

Neglecting the discretization of $K_{m,q,H}^t$ , $D$ and $\gamma_{q,H}$ , a series of time-split operators is defined by (17) [eq:vdfqs]). Once the diffusivities ( $K_{m,q,H}$ ) and the non-local transport terms ( $\gamma_{q,H}$ ) have been determined, the solution of (17) [eq:vdfqs]), proceeds in several steps.

update the bottom level values of $u$ , $v$ , $q$ and $s$ using the surface fluxes;
invert (17) and (19) for $u,v_{n+1}$ ;
compute $D$ and use to update the $s$ profile;
invert (17) [eq:vds]) and (20) for $s_{n+1}$ and $q_{n+1}$

Note that since all parameterizations in CAM4 return tendencies rather modified profiles, the actual quantities returned by the vertical diffusion are $\delta_n (u,v,s,q) / \delta t$ .

Equations (17) [eq:vdfqs]) constitute a set of four tridiagonal systems of the form

$-A^k \psi^{k+1}_{n+1} + B^k\psi^k_{n+1} -C^k\psi^{k-1}_{n+1} = \psi^k_{n\prime},$

where $\psi_{n\prime}$ indicates $u$ , $v$ , $q$ , or $s$ after updating from time $n$ values with the nonlocal and boundary fluxes. The super-diagonal ( $A^k$ ), diagonal ( $B^k$ ) and sub-diagonal ( $C^k$ ) elements of (21) are:

$\begin{aligned} A^k &=& \frac{1}{\delta^{k} p} \frac{\delta t}{\delta^{k+}p}\left(g^2\rho^2 K\right)_n^{k+},\\ B^k &=& 1 + A^k + C^k, \\ C^k &=& \frac{1}{\delta^{k} p} \frac{\delta t}{\delta^{k-}p}\left(g^2\rho^2 K\right)_n^{k-}.\end{aligned}$

The solution of (21) has the form

$\psi^k_{n+1} = E^k \psi^{k-1}_{n+1} + F^k,$

or,

$\psi^{k+1}_{n+1} = E^{k+1} \psi^{k}_{n+1} + F^{k+1}.$

Substituting (23) into (21) ,

$\psi^{k}_{n+1} = \frac{C^k}{B^k - A^k E^{k+1}} \psi^{k-1}_{n+1} + \frac{\psi^k_{n\prime} + A^k F^{k+1}}{B^k - A^k E^{k+1}}.$

Comparing (22) and (24) , we find

$E^k = \frac{C^k} {B^k - A^k E^{k+1}}, \quad L>k>1,$

$F^k = \frac{\psi^k_{n\prime} + A^k F^{k+1}}{B^k - A^k E^{k+1}}, \quad L>k>1$

The terms $E^k$ and $F^k$ can be determined upward from $k=L$ , using the boundary conditions

$E^{L+1} = F^{L+1} = A^L = 0.$

Finally, (24) can be solved downward for $\psi^{k}_{n+1}$ , using the boundary condition

$C^1 = 0 \Rightarrow E^1 = 0.$

CCM1-3 used the same solution method, but with the order of the solution reversed, which merely requires writing (23) for $\psi^{k-1}_{n+1}$ instead of $\psi^{k+1}_{n+1}$ . The order used here is particularly convenient because the turbulent diffusivities for heat and all constituents are the same but their molecular diffusivities are not. Since the terms in (25) and (26) are determined from the bottom upward, it is only necessary to recalculate $A^k$ , $C^k$ , $E^k$ and $1/({B^k - A^k E^{k+1}})$ for each constituent within the region where molecular diffusion is important.

7.6.3. Gravity Wave Drag¶

Vertically propagating gravity waves can be excited in the atmosphere where stably stratified air flows over an irregular lower boundary and by internal heating and shear. These waves are capable of transporting significant quantities of horizontal momentum between their source regions and regions where they are absorbed or dissipated. Previous GCM results have shown that the large-scale momentum sinks resulting from breaking gravity waves play an important role in determining the structure of the large-scale flow. CAM4 incorporates a parameterization for a spectrum of vertically propagating internal gravity waves based on the work of [Lin81], [Hol82], [GS85] and [McF87]. The parameterization solves separately for a general spectrum of monochromatic waves and for a single stationary wave generated by flow over orography, following [McF87]. The spectrum is omitted in the standard tropospheric version of CAM4, as in previous versions of the CCM. Here we describe the modified version of the gravity wave spectrum parameterization used in WACCM4.0.

7.6.3.1. Adiabatic inviscid formulation¶

Following (???), the continuous equations for the gravity wave parameterization are obtained from the two-dimensional hydrostatic momentum, continuity and thermodynamic equations in a vertical plane:

$\left( \frac{\partial}{\partial t} + u\frac{\partial}{\partial x}\right)u = -\frac{\partial\Phi}{\partial x}\, ,$

$\frac{\partial u}{\partial x} + \frac{\partial W}{\partial Z} = 0\, ,$

$\left( \frac{\partial}{\partial t} + u\frac{\partial}{\partial x}\right) \ \frac{\partial\Phi}{\partial Z} + N^2 w = 0\, .$

Where $N$ is the local Brunt-Väisällä frequency, and $W$ is the vertical velocity in log pressure height ( $Z$ ) coordinates. Eqs. (27) –(29) are linearized about a large scale background wind $\overline u$ , with perturbations $u^\prime,w^\prime$ , and combined to obtain:

$\left( \frac{\partial}{\partial t} + {\overline u}\frac{\partial}{\partial x}\right)^2 \frac{\partial^2 w^\prime}{\partial Z^2} + N^2 \frac{\partial^2 w^\prime}{\partial x^2} = 0\, .$

Solutions to (30) are assumed to be of the form:

$w^\prime = {\hat w} \, e^{ik(x-ct)} \, e^{Z/2H} \, ,$

where $H$ is the scale height, $k$ is the horizontal wavenumber and $c$ is the phase speed of the wave. Substituting (31) into (30) , one obtains:

$-k^2 (\overline u -c)^2 \left( \frac{\partial}{\partial Z} + \frac{1}{2H} \right)^2{\hat w} - k^2 N^2 {\hat w} = 0\, .$

Neglecting $\frac{1}{2H}$ compared to $\frac{\partial}{\partial Z}$ in eq:gw_what, one obtains the final form of the two dimensional wave equation:

$\frac{d^2 {\hat w}}{d Z^2} + \lambda^2 {\hat w} = 0\, ,$

with the coefficient defined as:

$\lambda= \frac{N}{(\overline u -c)}\, .$

The WKB solution of (33) is:

${\hat w} = A \lambda^{-1/2}\exp\left(i\int_0^Z\lambda dz^\prime\right)\, ,$

and the full solution, from (31) , is:

$w^\prime(Z,t) = A \lambda^{-1/2}\exp\left(i\int_0^Z\lambda dz^\prime\right) \ e^{ik(x-ct)}\ e^{Z/2H} \, .$

The constant $A$ is determined from the wave amplitude at the source ( $z=0$ ), The Reynolds stress associated with eq:eq:gw_final is:

$\tau(Z) = \tau(0) = \rho\overline{ u^\prime w^\prime} = -\frac{2}{k} |A|^2\rho_0{\rm sgn}(\lambda)\, ,$

and is conserved, while the momentum flux $\overline{ u^\prime w^\prime} = -(m/k)\ \overline{w^\prime w^\prime}$ grows exponentially with altitude as $\exp(Z/H)$ , per (36). We note that the vertical flux of wave energy is $c_{gz}\ E' = (U-c)\ \tau$ ((???)), where $c_{gz}$ is the vertical group velocity, so that deposition of wave momentum into the mean flow will be accompanied by a transfer of energy to the background state.

7.6.3.2. Saturation condition¶

The wave amplitude in (36) grows as $e^{Z/2H}$ until the wave becomes unstable to convective overturning, Kelvin-Helmholtz instability, or other nonlinear processes. At that point, the wave amplitude is assumed to be limited to the amplitude that would trigger the instability and the wave is “saturated”. The saturation condition used in CAM4 is from (???), based on a maximum Froude number ( $F_c$ ), or streamline slope.

$|\rho\overline{ u^\prime w^\prime}| \leq \tau^{*} = F_c^2\frac{k}{2}\rho \frac{|\overline{u}-c|^3}{N} \, ,$

where $\tau^*$ is the saturation stress and $F_c^2=0.5$ . In , $F_c^2=1$ and is omitted hereafter. Following (???), within a saturated region the momentum tendency can be determined analytically from the divergence of $\tau^*$ :

$\begin{aligned} \frac{\partial \overline u}{\partial t} &= -\frac{e}{\rho}\frac{\partial}{\partial Z} \rho\overline{ u^\prime w^\prime}\, , \nonumber \\ & \simeq -e \frac{k}{2} \frac{(\overline u-c)^3}{N} \frac{1}{\rho}\frac{\partial\rho}{\partial Z}\, , \nonumber \\ & \simeq -e \frac{k}{2} \frac{(\overline u-c)^3}{N H}, \end{aligned}$

where $e$ is an “efficiency” factor. For a background wave spectrum, $e$ represents the temporal and spatial intermittency in the wave sources. The analytic solution (39) is not used in WACCM4.0; it is shown here to illustrate how the acceleration due to breaking gravity waves depends on the intrinsic phase speed. In the model, the stress profile is computed at interfaces and differenced to get the specific force at layer midpoints.

7.6.3.3. Diffusive damping¶

In addition to breaking as a result of instability, vertically propagating waves can also be damped by molecular diffusion (both thermal and momentum) or by radiative cooling. Because the intrinsic periods of mesoscale gravity waves are short compared to IR relaxation time scales throughout the atmosphere, we ignore radiative damping. We take into account the molecular viscosity, $K_m^m$ , such that the stress profile is given by:

$\tau(Z) = \tau(Z_t) \exp\left(-\frac{2}{H}\int_0^Z\lambda_i dz^\prime\right) \, ,$

where $Z_t$ denotes the top of the region, below $Z$ , not affected by thermal dissipation or molecular diffusion. The imaginary part of the local vertical wavenumber, $\lambda_i$ is then:

$\lambda_i = \frac{N^3 \ K_m^m}{2 k (\overline u -c)^4 } \, .$

In WACCM4.0, ((40) – (41)) are only used within the domain where molecular diffusion is important (above $\sim 75$ km). At lower altitudes, molecular diffusion is negligible, $\lambda_i \rightarrow 0$ , and $\tau$ is conserved outside of saturation regions.

7.6.3.4. Transport due to dissipating waves¶

When the wave is dissipated, either through saturation or diffusive damping, there is a transfer of wave momentum and energy to the background state. In addition, a phase shift is introduced between the wave’s vertical velocity field and its temperature and constituent perturbations so that fluxes of heat and constituents are nonzero within the dissipation region. The nature of the phase shift and the resulting transport depends on the dissipation mechanism; in WACCM4.0, we assume that the dissipation can be represented by a linear damping on the potential temperature and constituent perturbations. For potential temperature, $\theta$ , this leads to:

$\left( {\partial \over \partial t} + \overline u {\partial \over \partial x} \right) \theta^\prime + w' {\partial \overline\theta \over \partial z} = -\delta \theta^\prime \, ,$

where $\delta$ is the dissipation rate implied by wave breaking, which depends on the wave’s group velocity, $c_{gz}$ (see [Gar01]).

$\delta = {c_{gz} \over 2H} = k \ {(\overline u - c)^2 \over 2H N} \, .$

Substitution of (43) into (42) then yields the eddy heat flux:

$\overline{w^\prime \theta^\prime} = -\left[ {\delta \ \overline{w^\prime w^\prime} \over k^2(\overline u - c)^2 + \delta^2} \right] {\partial \overline\theta \over \partial z} \, .$

Similar expressions can be derived for the flux of chemical constituents, with mixing ratio substituted in place of potential temperature in ([eq:gw:sub:kzz]). We note that these wave fluxes are always downgradient and that, for convenience of solution, they may be represented as vertical diffusion, with coefficient $K_{zz}$ equal to the term in brackets in (44) , but they do not represent turbulent diffusive fluxes but rather eddy fluxes. Any additional turbulent fluxes due to wave breaking are ignored. To take into account the effect of localization of turbulence (e.g., [FD85][McI89]) (44) is multiplied times an inverse Prandtl number, ${Pr}^{-1}$ ; in WACCM4.0 we use ${Pr}^{-1}=0.25$ .

7.6.3.5. Heating due to wave dissipation¶

The vertical flux of wave energy density, $E'$ , is related to the stress according to:

$c_{gz} \ E' = (\overline u-c) \ \tau \ ,$

where $c_{gz}$ is the vertical group velocity [Andrews et al., 1987]. Therefore, the stress divergence $\partial \tau / \partial Z$ that accompanies wave breaking implies a loss of wave energy. The rate of dissipation of wave energy density is:

${\partial E' \over \partial t} \simeq (\overline u-c) {1 \over c_{gz}}{\partial \tau \over \partial t} = (\overline u-c) {\partial \tau \over \partial Z} \ .$

For a saturated wave, the stress divergence is given by (39) , so that:

${\partial E' \over \partial t} = (\overline u - c) \ {\partial \, \tau^* \over \partial Z} = - e \cdot \rho \ {k \, (U-c)^4 \over 2 N H} \ .$

This energy loss by the wave represents a heat source for the background state, as does the change in the background kinetic energy density implied by wave drag on the background flow:

${\partial \overline K \over \partial t} \equiv {\rho \over 2} {\partial \overline u^2 \over \partial t} = \overline u \ {\partial \, \tau^* \over \partial Z} = -e \cdot \rho \ {k \, \overline u \, (\overline u-c)^3 \over 2 NH} \ ,$

which follows directly from (39) . The background heating rate, in K sec $^{-1}$ , is then:

$Q_{gw} = -{1 \over \rho\, c_p} \left[{\partial \overline K \over \partial t} + {\partial E' \over \partial t} \right].$

Using $(\ref{eq:gw_eprime})-(\ref{eq:gw_kbar})$ , this heating rate may be expressed as:

$Q_{gw} = {1 \over \rho\, c_p} \ c \ {\partial \, \tau^* \over \partial Z} = {1 \over c_p} \left[ \ e \cdot {k \, c\,(c-\overline u)^3 \over 2 N H} \right] ,$

where $c_p$ is the specific heat at constant pressure. In WACCM4.0, $Q_{gw}$ is calculated for each component of the gravity wave spectrum using the first equality in (47) , i.e., the product of the phase velocity times the stress divergence.

7.6.3.6. Orographic source function¶

For orographically generated waves, the source is taken from (???):

$\tau_g = |\rho\overline{ u^\prime w^\prime}|_0 = \frac{k}{2} h_0^2 \rho_0 N_0 \overline u_0\, ,$

where $h_0$ is the streamline displacement at the source level, and $\rho_0$ , $N_0$ , and $\overline u_0$ are also defined at the source level. For orographic waves, the subgrid-scale standard deviation of the orography $\sigma$ is used to estimate the average mountain height, determining the typical streamline displacement. An upper bound is used on the displacement (equivalent to defining a “separation streamline”) which corresponds to requiring that the wave not be supersaturated at the source level:

$h_0=\min(2\sigma, \frac{\overline u_0}{N_0})\, .$

The source level quantities $\rho_0$ , $N_0$ , and $\overline u_0$ are defined by vertical averages over the source region, taken to be $2\sigma$ , the depth to which the average mountain penetrates into the domain:

$\psi_0 = \int_0^{2\sigma} \psi\rho dz, \qquad \psi \in \{\rho,N, u, v\} \, .$

The source level wind vector $(u_0,v_0)$ determines the orientation of the coordinate system in (27) –(29) and the magnitude of the source wind $\overline u_0$ .

7.6.3.7. Non-orographic source functions¶

The source spectrum for non-orographic gravity waves is no longer assumed to be a specified function of location and season, as was the case with the earlier version of the model described by (???). Instead, gravity waves are launched according to trigger functions that depend on the atmospheric state computed in WACCM4 at any given time and location, as discussed by (???). Two trigger functions are used: convective heat release (which is a calculated model field) and a “frontogenesis function”, (???), which diagnoses regions of strong wind field deformation and temperature gradient using the horizontal wind components and potential temperature field calculated by the model.

In the case of convective excitation, the method of (???) is used to determine a phase speed spectrum based upon the properties of the convective heating field. A spectrum is launched whenever the deep convection parameterization in WACCM4 is active, and the vertical profile of the convective heating, together with the mean wind field in the heating region, are used to determine the phase speed spectrum of the momentum flux. Convectively generated waves are launched at the top of the convective region (which varies according to the depth of the convective heating calculated in the model).

Waves excited by frontal systems are launched whenever the frontogenesis trigger function exceeds a critical value (see (???)). The waves are launched from a constant source level, which is specified to be 600 mb. The momentum flux phase speed spectrum is given by a Gaussian function in phase speed:

$\tau_s (c) = \tau_b \exp\left[-\left( \frac{c - V_s}{c_w} \right)^2\right],$

centered on the source wind, $V_s = |{\bf V}_s|$ , with width $c_w=30$ m/s. A range of phase speeds with specified width and resolution is used:

$c \in V_s + [\pm d_c, \pm 2d_c, ... \pm c_{max}]\, ,$

with $d_c = 2.5$ m s $^{-1}$ and $c_{max} = 80$ m s $^{-1}$ , giving 64 phase speeds. Note that $c=V_s$ is retained in the code for simplicity, but has a critical level at the source and, therefore, $\tau_s(c=V_s)=0$ . Note also that $\tau_b$ is a tunable parameter; in practice this is set such that the height of the polar mesopause, which is very sensitive to gravity wave driving, is consistent with observations. In WACCM4, $\tau_b$ = 1.5 x 10 $^{-3}$ Pa.

Above the source region, the saturation condition is enforced separately for each phase speed, $c_i$ , in the momentum flux spectrum:

$\tau(c_i) \leq \tau_i^{*} =F_c^2 \frac{k}{2}\rho \frac{|\overline u-c_i|^3}{N} .$

7.6.3.8. Numerical approximations¶

The gravity wave drag parameterization is applied immediately after the nonlinear vertical diffusion. The interface Brunt-Väisällä frequency is

$\left(N^{k+}\right)^2 = \frac{g^2}{T^{k+}} \left( \frac{1}{c_p} - \rho^{k+} \frac{\delta^{k+} T}{\delta^{k+} p} \right) \, ,$

Where the interface density is:

$\rho^{k+} = \frac{R T^{k+}}{p^{k+}}.$

The midpoint Brunt-Väisällä frequencies are $N^k=(N^{k+}+N^{k-})/2$ .

The level for the orographic source is an interface determined from an estimate of the vertical penetration of the subgrid mountains within the grid box. The subgrid scale standard deviation of the orography, $\sigma_h$ , gives the variation of the mountains about the mean elevation, which defines the Earth’s surface in the model. Therefore the source level is defined as the interface, $k_s-1/2$ , for which $z^{k_s+} < 2\sigma_h < z^{k_s-}$ , where the interface heights are defined from the midpoint heights by $z^{k+} = \sqrt(z^{k} z^{k+1})$ .

The source level wind vector, density and Brunt-Väisällä frequency are determined by vertical integration over the region from the surface to interface $k_s+1/2$ :

$\psi_0 = \sum_{k=k_s}^K \psi^k \delta^k p\, , \qquad \psi \in \{\rho,N, u, v\} \, .$

The source level background wind is $\overline u_0 =\sqrt(u_0^2 + v_0^2)$ , the unit vector for the source wind is

$(x_0, y_0) = (u_0, v_0) /\overline u_0 \, ,$

and the projection of the midpoint winds onto the source wind is

$\overline u^k = u^k x_0 + v^k y_0 \, .$

Assuming that $\overline u_0 > 2$ m s $^{-1}$ and $2\sigma^h > 10$ m, then the orographic source term, $\tau_g$ is given by (48) and (49) , with $F_c^2$ =1 and $k = 2\pi/10^5$ m $^{-1}$ . Although the code contains a provision for a linear stress profile within a “low level deposition region”, this part of the code is not used in the standard model.

The stress profiles are determined by scanning up from the bottom of the model to the top. The stress at the source level is determined by (48) . The saturation stress, $\tau^*_\ell$ at each interface is determined by (53) , and $\tau^*_\ell=0$ if a critical level is passed. A critical level is contained within a layer if $(\overline u^{k+} -c_\ell) / (\overline u^{k-} -c_\ell) < 0$ .

Within the molecular diffusion domain, the imaginary part of the vertical wavenumber is given by (41) . The interface stress is then determined from the stress on the interface below by:

$\tau^{k-} = \min \left[ \left(\tau^*\right)^{k-}, \tau^{k+}\exp\left( -2 \lambda_i \frac{R}{g} T^k \delta^k\ln p \right) \right] \, .$

Below the molecular diffusion domain, the exponential term in (58) is omitted.

Once the complete stress profile has been obtained, the forcing of the background wind is determined by differentiating the profile during a downward scan:

$\frac{\partial \overline u^k_\ell}{\partial t} = g \frac{\delta^k\tau_\ell}{\delta^k p} < \left( \frac{\partial \overline u^k_\ell}{\partial t}\right)^{\rm max} \, .$

$\left( \frac{\partial \overline u^k_\ell}{\partial t}\right)^{\rm max} = \min\left[ \frac{|c_\ell - \overline u^k_\ell|}{2 \delta t}\, , 500 {\rm ~m ~s^{-1}~ day^{-1}} \right] \, .$

The first bound on the forcing comes from requiring that the forcing not be large enough to push the wind more than half way towards a critical level within a time step and takes the place of an implicit solution. This bound is present for numerical stability, it comes into play when the time step is too large for the forcing. It is not feasible to change the time step, or to write an implicit solver, so an a priori bound is used instead. The second bound is used to constrain the forcing to lie within a physically plausible range (although the value used is extremely large) and is rarely invoked.

When any of the bounds in (59) are invoked, conservation of stress is violated. In this case, stress conservation is ensured by decreasing the stress on the lower interface to match the actual stress divergence in the layer:

$\tau^{k+}_\ell = \tau^{k-}_\ell + \frac{\partial \overline u^k}{\partial t} \frac{\delta^k p}{g} \, .$

This has the effect of pushing some of the stress divergence into the layer below, a reasonable choice since the waves are propagating up from below.

Finally, the vector momentum forcing by the gravity waves is determined by projecting the background wind forcing with the unit vectors of the source wind:

$\frac{\partial {\bf V}^k}{\partial t} = (x_0, y_0) \times E \sum_\ell \frac{\partial \overline u^k_\ell}{\partial t} \, .$

In addition, the frictional heating implied by the momentum tendencies, $\frac{1}{c_p} {\bf V}^k \cdot {\partial {\bf V}^k / \partial t}$ , is added to the thermodynamic equation. This is the correct heating for orographic ( $c_\ell=0$ ) waves, but not for waves with $c_\ell\ne 0$ , since it does not account for the wave energy flux. This flux is accounted for in some middle and upper atmosphere versions of CAM4, but also requires accounting for the energy flux at the source.

7.6.4. Turbulent Mountain Stress¶

An important difference between WACCM4 and earlier versions is the addition of surface stress due to unresolved orography. A numerical model can compute explicitly only surface stresses due to resolved orography. At the standard 1.9 $^\circ$ x 2.5 $^\circ$ (longitude x latitude) resolution used by WACCM4 only the gross outlines of major mountain ranges are resolved. To address this problem, unresolved orography is parameterized as turbulent surface drag, using the concept of effective roughness length developed by (???). Fiedler and Panofsky defined the roughness length for heterogeneous terrain as the roughness length that homogenous terrain would have to give the correct surface stress over a given area. The concept of effective roughness has been used in several Numerical Weather Prediction models (e.g., (???); (???)).

In WACCM4 the effective roughness stress is expressed as:

$\tau = \rho \, C_d \, |{\bf V}| {\bf V} \, ,$

where $\rho$ is the density and $C_d$ is a turbulent drag coefficient,

$C_d = \frac{f(R_i)\, k ^2}{\ln^2\left[\frac{z+z_0}{z_0}\right]} \, ,$

$k$ is von Kármán’s constant; $z$ is the height above the surface; $z_0$ is an effective roughness length, defined in terms of the standard deviation of unresolved orography; and $f(R_i)$ is a function of the Richardson number (see (???) for details).

The stress calculated by (63) is used the model’s nonlocal PBL scheme to evaluate the PBL height and nonlocal transport, per Eqs. (3.10)Ð(3.12) of (???). This calculation is carried out only over land, and only in grid cells where the height of topography above sea level, $z$ , is nonzero.

7.6.5. QBO Forcing¶

WACCM4 has several options for forcing a quasi-biennial oscillation (QBO) by applying a momentum forcing in the tropical stratosphere. The parameterization relaxes the simulated winds to a specified wind field that is either fixed or varies with time. The parameterization can also be turned off completely. The namelist variables and input files can be selected to choose one of the following options:

Idealized QBO East winds, used for perpetual fixed-phase of the QBO, as described by (???).
Idealized QBO West winds, as above but for the west phase.
Repeating idealized 28-month QBO, also described by (???).
QBO for the years 1953-2004 based on the climatology of Giorgetta [see: http://www.pa.op.dlr.de/CCMVal/Forcings/qbo_data_ccmval/u_profile_195301-200412.html, 2004].
QBO with a 51-year repetition, based on the 1953-2004 climatology of Giorgetta, which can be used for any calendar year, past or future.

The relaxation of the zonal wind is based on (???) and is described in (???). The input winds are specified at the equator and the parameterization extends latitudinally from 22 $^{\circ}$ N to 22 $^{\circ}$ S, as a Gaussian function with a half width of 10 $^{\circ}$ centered at the equator. Full vertical relaxation extends from 86 to 4 hPa with a time constant of 10 days. One model level below and above this altitude range, the relaxation is half as strong and is zero for all other levels. This procedure constrains the equatorial winds to more realistic values while allowing resolved and parameterized waves to continue to propagate.

The fixed or idealized QBO winds (first 3 options) can be applied for any calendar period. The observed input (Giorgetta climatology) can be used only for the model years 1953-2004. The winds in the final option were determined from the Giorgetta climatology for 1954-2004 via filtered spectral decomposition of that climatology. This gives a set of Fourier coefficients that can be expanded for any day and year. The expanded wind fields match the climatology during the years 1954-2004.

7.6.6. Radiation¶

The radiation parameterizations in CAM4 are quite accurate up to $\sim 65$ km, but deteriorate rapidly above that altitude. Because 65 km is near a local minimum in both shortwave heating and longwave cooling, it is a particularly convenient height to merge the heating rates from parameterizations for the lower and upper atmosphere. Therefore, we retain the CAM4 parameterizations below $\sim 65$ km and use new parameterizations above.

The merged shortwave and longwave radiative heatings are determined from

$Q = w_1 \, Q_{CAM3} + w_2 \, Q_{MLT},$

where $w_1(z^*<z_b^*) = 1$ , $w_2(z^*>z_t^*) = 1$ and $z^* = \log(10^5/p)$ is the pressure scale height. The CAM4 radiation parameterizations are used below $z_b^*$ and the MLT parameterizations are used above $z_t^*$ . For $z_b^* < z < z_t^*$ , $w_2 = 1 - w_1$ and

$w_1 = 1 - \tanh\left( \frac{z^* - z_b^*}{z_w*}\right),$

where $z_w*$ is the transition width.

The merging was developed and tested separately for shortwave and longwave radiation and the constants are slightly different. For longwave radiation, the constants are $z_b^*=8.57$ , $z_t^*=10$ and $z_w^*=0.71$ . For shortwave radiation, the constants are $z_b^*=9$ , $z_t^*=10$ and $z_w^*=0.75$ . These constants give smooth heating profiles. Note that a typical atmospheric scale height of $H=7$ km places the transition zones between 60 and 70 km.

7.6.6.1. Longwave radiation¶

WACCM4.0 retains the longwave (LW) formulation used in CAM4 (???). However, in the MLT longwave radiation uses the parameterization of (???) for $\rm CO_2$ and $\rm O_3$ cooling and the parameterization of (???) for $\rm NO$ cooling at 5.3 $\mu$ m. As noted above, the LW heating/cooling rates produced by these parameterizations are merged smoothly at 65 km with those produced by the standard CAM4 LW code, as recently revised by (???). In the interactive chemistry case all of the gases (O, $\rm O_2$ , $\rm O_3$ , $\rm N_2$ , NO, and $\rm CO_2$ ) that are required by these parameterizations, are predicted within WACCM4.0.

7.6.6.2. Shortwave radiation¶

WACCM4.0 uses a combination of solar parameterizations to specify spectral irradiances over two spectral intervals. The first spectral interval covers soft x-ray and extreme ultraviolet irradiances (wavelengths between 0.05 nm to Lyman- $\alpha$ (121.6 nm)) and is calculated using the parameterization of (???). The parameterizations take as input the 10.7 cm solar radio flux ( $f10.7$ ) and its 81-day average ( $f10.7a$ ). Daily values of $f10.7$ are obtained from NOAA’s Space Environment Center (www.sec.noaa.gov).

The irradiance of the $j$ th spectral interval is:

$F_j = F_j^0 * \left\{ 1 + R_j*\left[\frac{(f10.7 + f10.7a)}{2}-F_{min}\right] \right\}$

where $F_{min}$ = 80. $F_j^0$ and $R_j$ are taken from Table A1 of (???).

Fluxes for the second interval between Lyman- $\alpha$ (121.6 nm) and 100 $\mu$ m. are specified using an empirical model of the wavelength-dependent sunspot and facular influences (???; Wang, Lean, and Sheeley 2005). Spectral resolution is 1 nm between 121.6 nm and 750nm, 5 nm between 750nm and 5 $\mu$ m, 10 nm between 5 $\mu$ m and 10 $\mu$ m, and 50 nm between 10 $\mu$ m and 100 $\mu$ m.

In the troposphere, stratosphere and lower mesosphere (z $<$ 65km) WACCM4.0 retains the CAM4 shortwave heating (200 nm to 4.55 $\mu$ m) which is calculated from the net shortwave spectral flux into each layer (???). The solar spectrum for the CAM4 heating calculation is divided into 19 intervals (???). The heating in these intervals must be adjusted to match the irradiances calculated for the upper part of the model, and those used in the photolysis calculations. This is achieved by applying a scaling ( $S_j$ ) to the solar heating in the $j$ th CAM4 spectral interval using the spectrum from (???) and Wang, Lean, and Sheeley (2005):

$S_j = \frac{F_{j}}{F_{j}^{ref}},$

where $F_j$ is the spectral irradiance (W/m:math:^2/nm) integrated over the $j$ th band, and $F_{j}^{ref}$ is the same integral taken over a reference spectrum calculated from annual mean fluxes over a 3-solar-cycle period from XX to YY.

In the MLT region, shortwave heating is the sum of the heating due to absorption of photons and subsequent exothermic chemical reactions that are initiated by photolysis. The majority of energy deposited by an absorbed photon goes into breaking molecular bonds, rather than into translational energy of the absorbing molecule (heat). Chemical heating results when constituents react to form products of lower total chemical potential energy. This heating can take place months after the original photon absorption and thousands of kilometers away. Heating rates range from 1 K/day near 75 km to 100-300 K/day near the top of the model domain. It is clear that quenching of $O(^1D)$ is a large source of heating throughout the MLT. Above 100 km ion reactions and reactions involving atomic nitrogen are significant sources of heat, while below that level O $_X$ (= O + O $_3$ ) and HO $_X$ (= H + OH + HO $_2$ ) reactions are the dominant producers of chemical heating.

Heating within the MLT from the absorption of radiation that is directly thermalized is calculated over the wavelength range of 0.05 nm to 350 nm. For wavelengths less than Lyman- $\alpha$ , it is assumed that 5% of the energy of each absorbed photon is directly thermalized:

$Q_{EUV} = (\rho c_p)^{-1} \sum_{k} n_k \sum_{j} \epsilon J_{k}(\lambda_j) \frac{hc}{\lambda_j},$

where $\epsilon$ = 0.05. Here $\rho$ is mass density, $c_p$ is the specific heat of dry air, $n$ is the number density of the absorbing species, and $J$ is the photolysis/photoionization rate. The total heating is the sum of $k$ photolysis reactions and $j$ wavelengths intervals. At these wavelengths absorption of a photon typically leads to photoionization, with the resulting photoelectron having sufficient energy to ionize further molecules. Calculation of $J_{ij}$ and ionization rates from photoelectrons is calculated based on the parameterization of (???). In a similar manner, the heating rate within the aurora ( $Q_{AUR}$ ) is calculated as the product of the total ionization rate, 35 eV per ion pair, and the same heating efficiency of 5%.

Between Lyman- $\alpha$ and 350 nm the energy required to break molecular bonds is explicitly accounted for. The heating rate is thus defined as:

$Q_{UV} = (\rho c_p)^{-1} \sum_{k} n_k \sum_{j} J_{k}(\lambda_j) \{ \frac{hc}{\lambda_j} -BDE_k \},$

where $BDE$ is the bond dissociation energy.

In addition to these sources of heat, WACCM4.0 calculates heating by absorption in the near-infrared by CO $_2$ (between 1.05 to 4.3 $\mu$ m), which has its largest contribution near 70km and can exceed 1 K/day (???). Heating from this process is calculated using the parameterization of (???). Finally, the heating produced by collisions of electrons and neutrals (Joule heating) is also calculated using the predicted ion and electron concentrations. This is described in section [ion:sub:drag]. Local heating rates from joule heating can be very large in the auroral regions, reaching over 10 $^3$ K/day in the upper levels of the model.

Airglow, radiation produced when excited atoms or molecules spontaneously emit, is accounted for in WACCM4.0 for emissions of O $_2(^1\Delta )$ , O $_2(^1\Sigma )$ , and vibrationally excited OH. Airglow from the excited molecular oxygen species are handled explicitly; radiative lifetimes for O $_2(^1\Delta )$ and O $_2(^1\Sigma )$ are 2.58 $\times$ 10 $^{-4}$ s $^{-1}$ and 0.085 s $^{-1}$ respectively. However, modeling of the many possible vibrational transitions of OH is impractical in a model as large as WACCM4.0. Energy losses from the emission of vibrationally excited OH are therefore accounted for by applying an efficiency factor to the exothermicity of the reaction that produces vibrationally excited OH; the reaction of hydrogen and ozone. In other words, the reaction H + O $_3$ produces ground state OH only, but the chemical heating from the reaction has been reduced to take into consideration that some of the chemical potential energy has been lost in airglow. This approach is the same one used by (???) and we use their recommended efficiency factor of 60%. Any energy lost through airglow is assumed to be lost to space, and so represents an energy pathway that does not generate heat.

7.6.6.3. Volcanic Heating¶

The sulfate aerosol heating is a function of a prescribed aerosol distribution varying in space and time that has a size distribution similar to that seen after a volcanic eruption (???). The H $_{\rm 2}$ SO $_{\rm 4}$ mass distribution is calculated from the prescribed sulfate surface area density (SAD) assuming a lognormal size distribution, number of particles per cm-3, and distribution width (see section 3.6.2). The H2SO4 mass distribution is then passed to the radiative transfer code (CAMRT), which in turn calculates heating and cooling rates.

7.6.7. chemistry¶

7.6.7.1. Chemical Mechanism (Neutral Species)¶

includes a detailed neutral chemistry model for the middle atmosphere based on the Model for Ozone and Related Chemical Tracers, Version 3 (???). The mechanism represents chemical and physical processes in the troposphere through the lower thermosphere. The species included within this mechanism are contained within the O $_{\rm X}$ , NO $_{\rm X}$ , HO $_{\rm X}$ , ClO $_{\rm X}$ , and BrO $_{\rm X}$ chemical families, along with CH $_4$ and its degradation products. This mechanism contains 52 neutral species, one invariant (N:math:_2), 127 neutral gas-phase reactions, 48 neutral photolytic reactions, and 17 heterogeneous reactions on three aerosol types (see below). Lists of the chemical species are given in Table 1. The first column lists the symbolic name (as used in the mechanism); the second column lists the species atomic composition; the third column designates which numerical solution approach is used (i.e., explicit or implicit); the fourth column lists any deposition processes (wet or dry) for that species; and the fifth column indicates whether the surface (or upper) boundary condition is fixed vmr or flux, or if a species has an in-situ flux (from lightning or aircraft emissions).

The gas-phase reactions included in the middle atmosphere chemical mechanism are listed in Table 2. In most all cases the chemical rate constants are taken from JPL06-2 (???). Exceptions to this condition are described in the comment section for any given reaction.

Heterogeneous reactions on four different aerosols types are also represented in the chemical mechanism (see Table 3): 1) liquid binary sulfate (LBS); 2) Supercooled ternary solution (STS); 3) Nitric acid trihydrate (NAT); and 4) water-ice. There are 17 reactions, six reactions on liquid sulfate aerosols (LBS or STS), five reactions on solid NAT aerosols, and six reactions on solid water-ice aerosols. The rate constants for these 17 heterogeneous reactions can be divided up into two types: 1) first order; and 2) pseudo second order. For first order hydrolysis reactions (Table 3, reactions 1-3, 7-8, 11, and 12-14), the heterogeneous rate constant is derived in the following manner:

$k=\frac{1}{4}V\cdot SAD \cdot \gamma$

Where V = mean velocity; SAD = surface area density of LBS, STS, NAT, or water-ice, and $\gamma$ = reaction probability for each reaction. The units for this rate constant are s $^{-1}$ . Here the H $_2$ O abundance is in excess and assumed not change relative to the other reactant trace constituents. The mean velocity is dependent on the molecular weight of the non-H $_2$ O reactant (i.e., N $_2$ O $_5$ , ClONO $_2$ , or BrONO $_2$ ). The SAD for each aerosol type is described in section 7. The reaction probability is dependent on both composition and temperature for sulfate aerosol (see JPL06-2). The reaction probability is a fixed quantity for NAT and water-ice aerosols and is listed in Table 3. Multiplying the rate constant times the concentration gives a loss rate in units of molecules cm $^{-3}$ sec $^{-1}$ for the reactants and is used in the implicit solution approach. The non-hydrolysis reaction (Table 3, reactions 4-6, 9-10, and 15-17) are second order reactions. Here, the first order rate constant (equation 6) is divided by the HCl concentration, giving it the typical bimolecular rate constant unit value of cm $^3$ molecule $^{-1}$ sec $^{-1}$ . This approach assumes that all the HCl is in the aerosol particle.

7.6.7.2. Stratospheric Aerosols¶

Heterogeneous processes on liquid sulfate aerosols and solid polar stratospheric clouds (Type 1a, 1b, and 2) are included following the approach of (???). This approach assumes that the condensed phase mass follows a lognormal size distribution taken from (???),

$N(r) = \frac{N_0}{r \sigma \sqrt{2\pi}}\exp\left[\frac{-\ln (r / r_0)^2} {2 \sigma ^2}\right]$

where $N$ is the aerosol number density (particles cm $^{-3}$ ); $r$ and $r_0$ are the particle radius and median radius respectively; and $\sigma$ is the standard deviation of the lognormal distribution. $N_0$ and $r_0$ are supplied for each aerosol type. The aerosol surface area density (SAD) is the second moment of this distribution.

At model temperatures (T:math:_{rm model}) greater than 200 K, liquid binary sulfate (LBS) is the only aerosol present. The surface area density (SAD) for LBS is derived from observations from SAGE, SAGE-II and SAMS (???) as updated by Considine (???). As the model atmosphere cools, the LBS aerosol swells, taking up both HNO $_3$ and H $_2$ O to give STS aerosol. The Aerosol Physical Chemistry Model (ACPM) is used to derive STS composition (???). The STS aerosol median radius and surface area density is derived following the approach of (???). The width of the STS size distribution ( $\sigma=1.6$ ) and number density (10 particles cm $^{-3}$ ) are prescribed according to measurements from (???). The STS aerosol median radius can swell from approximately 0.1 $\mu$ m to approximately 0.5 $\mu$ m. There is no aerosol settling assumed for this type of aerosol. The median radius is used in derivation of sulfate aerosol reaction probability coefficients. Both the LBS and STS surface area densities are used for the calculation of the rate constants as listed in Table 3; reactions (1)-(6).

Solid nitric acid containing aerosol formation is allowed when the model temperature reaches a prescribed super saturation ratio of HNO $_3$ over NAT [Hansen and Mauersberger, 1988]. This ratio is set to 10 in (???). There are three methods available to handle the HNO $_3$ uptake on solid aerosol. The first method directly follows (???; ???). Here, after the supersaturation ratio assumption is met, the available condensed phase HNO $_3$ is assumed to reside in the solid NAT aerosol. The derivation of the NAT median radius and surface area density follows the same approach as the STS aerosol, by assuming: a lognormal size distribution, a width of a distribution ( $\sigma=1.6$ ; (???)), and a number density (0.01 particles cm $^{-3}$ ; (???)). The NAT radius settles with a value of $r_0$ ranging between 2 and 5 $\mu$ m; this value depends on the model temperature and subsequent amount of condensed phase HNO $_3$ formed. This NAT median radius $r_0$ is also used to derive the terminal velocity for settling of NAT (section 8) and the eventual irreversible denitrification. The NAT surface area density is used to calculate the rate constants for heterogeneous reactions 7-11 (Table 3). Since the available HNO $_3$ is included inside the NAT aerosol, there is no STS aerosol present. However, there are still heterogeneous reactions occurring on the surface of LBS aerosols.

If the calculated atmospheric temperature, $T$ , becomes less than or equal to the saturation temperature (T:math:_{sat}) for water vapor over ice (e.g., (???)), water-ice aerosols can form. In the condensed phase H $_2$ O is derived in the prognotic water routines of CAM and passed into the chemistry module. Using this condensed phase H $_2$ O, the median radius and the surface area density for water-ice are again derived following the approach of (???). The water-ice median radius and surface area density assumes a lognormal size distribution, a width of a distribution = 1.6 (???), and a number density of 0.001 particles cm $^{-3}$ (???). The value of $r_0$ is typically 10 $\mu$ m. The water-ice surface area density is used for the calculation of the rate constants for reactions 12-17 (Table 3).

7.6.7.3. Sedimentation of Stratospheric Aerosols¶

The sedimentation of HNO $_3$ in stratospheric aerosols follows the approach described in (???). The following equation is used to derive the flux ( $F$ ) of HNO $_3$ , as NAT aerosol, across model levels in units of molecules cm $^{-2}$ sec $^{-1}$ .

$F_i = V_i \cdot C_i \, \exp(8 \ln^2\sigma_i),$

where $i = 1$ for NAT; $V_i$ is the terminal velocity of the aerosol particles (cm s $^{-1}$ ); $C$ is the condensed-phase concentration of HNO $_3$ (molecules cm $^{-3}$ ); $\sigma$ is the width of the lognormal size distribution for NAT (see discussion above). The terminal velocity is dependent on the given aerosol: 1) mass density; 2) median radius; 3) shape; 4) dynamic viscosity; and 5) Cunningham correction factor for spherical particles (see (???) and (???) for the theory behind the derivation of terminal velocity). For each aerosol type the terminal velocity could be calculated, however, in this quantity is only derived for NAT. Settling of HNO $_3$ contain in STS is not derived based on the assumption that the median radius is too small to cause any significant denitrification and settling of condensed phase H $_2$ O is handled in the CAM4 prognostic water routines.

7.6.7.4. Ion Chemistry¶

includes a six constituent ion chemistry model (O:math:^+, O $_2^+$ , N $^+$ , N $_2^+$ , NO $^+$ , and electrons) that represents the the E-region ionosphere. The global mean ion and electron distributions simulated by for solar minimum conditions are shown in Figure [ionfig], which clearly shows that the dominant ions in this region are NO $^+$ and O $_2^+$ . Ion-neutral and recombination reactions included in are listed in Table [ionrxntab]. The reaction rate constants for these reactions are taken from (???).

Ionization sources include not only the aforementioned absorption of extreme ultraviolet and soft x-ray photons, and photoelectron impact, but also energetic particles precipitation in the auroral regions. The latter is calculated by a parameterization based on code from the NCAR TIME-GCM model (???) that rapidly calculates ion-pair production rates, including production in the polar cusp and polar cap. The parameterization takes as input hemispheric power (HP), the estimated power in gigawatts deposited in the polar regions by energetic particles.

Currently uses a parameterization of HP (in GW) based on an empirical relationships between HP and the $K_p$ planetary geomagnetic index. For $K_p \leq 7$ , uses the relationship obtained by (???) from TIMED/GUVI observations:

${\rm HP} = 16.82*K_p*\exp{(0.32)}-4.86$

For $K_p > 7$ , linearly interpolates HP, assuming HP equals to 300 when $K_p$ is 9, based on NOAA satellite measurements:

${\rm HP} = 153.13 + \frac{K_p - 7}{9-7}*(300 - 153.13)$

$K_p$ is also available from NOAA’s Space Environment Center and covers the period from 1933 to the present, making it ideal for long-term retrospective simulations.

Global mean distribution of charged constituents during July solar minimum conditions.

a) Global distribution of ionization rates at 7.3 $\times 10^{-5}$ hPa, July 1, UT0100 HRS. Contour interval is 2 $\times 10^{3}$ cm $^{-3}$ s $^{-1}$ . b) Simultaneous global mean ionization rates (cm:math:^{-3} s $^{-1}$ ) versus pressure.

Total ionization rates at 110km during July for solar maximum conditions are shown in Figure [qsumfig]a. The broad region of ionization centered in the tropics is a result of EUV ionization, and has a peak value of almost 10 $^3$ at 22 $^\circ$ N. Ionization rates from particle precipitation can exceed this rate by 40% but are limited to the high-latitudes, as can been seen by the two bands that are approximately aligned around the magnetic poles. The global mean ionization rate (Figure [qsumfig]b)

An important aspect of including ionization processes (both in the aurora and by energetic photons and photoelectrons), is that it leads to a more accurate representation of thermospheric nitric oxide. Not only does nitric oxide play an important role in the energy balance of the lower thermosphere through emission at 5.3 $\mu$ m, it might also be transported to the upper stratosphere, where it can affect ozone concentrations. Nitric oxide is produced through quenching of N $(^2D)$ :

$N(^2D) + O_2 \rightarrow NO + O(^1D) + 1.84eV$

N $(^2D)$ is produced either via recombination of NO $^+$ (see Table [ionrxntab]) or directly by ionization of molecular nitrogen. The branching ratio between N $(^2D)$ and ground-state atomic nitrogen for the photoionization process is critical in determining the effectiveness of NO production. If ground-state atomic nitrogen is produced then it can react with NO to produce molecular nitrogen and effectively remove to members of the NOx family. In 60% of the atomic nitrogen produced is in the excited state, which implies absorption of EUV results in a net source of NO. Also shown are maxima at high latitudes due to auroral ionization. reproduces many of the features of the Nitric Oxide Empirical Model (NOEM) distribution (???), which is based on data from the Student Nitric Oxide Explorer satellite (???) In particular, larger NO in the winter hemisphere (a result of less photolytic loss), and a more localized NO maximum in the Northern Hemisphere (related to the lesser offset of geographic and magnetic poles, and so less spread when viewed as a geographic zonal mean).

no.	Symbolic Name	Chemical Formula	Numerics	Deposition	Boundary Condition
1	O	O( $^3$ P)	Implicit		ubvmr
2	O1D	O( $^1$ D)	Implicit
3	O3	O $_3$	Implicit	dry
4	O2	O $_2$	Implicit		ubvmr
5	O2_1S	O $_2(^1\Sigma)$	Implicit
6	O2_1D	O $_2(^1\Delta)$	Implicit
7	H	H	Implicit		ubvmr
8	OH	OH	Implicit
9	HO2	HO $_2$	Implicit
10	H2	H $_2$	Implicit		vmr, ubvmr
11	H2O2	H $_2$ O $_2$	Implicit	dry, wet
12	N	N	Implicit		ubvmr
13	N2D	N( $^2$ D)	Implicit		from TIME-GCM
14	N2	N $_2$	Invariant
15	NO	NO	Implicit		flux, ubvmr,
					lflux, airflux
16	NO2	NO $_2$	Implicit	dry
17	NO3	NO $_3$	Implicit
18	N2O5	N $_2$ O $_5$	Implicit
19	HNO3	HNO $_3$	Implicit	dry, wet
20	HO2NO2	HO $_2$ NO $_2$	Implicit	dry, wet
21	CL	Cl	Implicit
22	CLO	ClO	Implicit
23	CL2	Cl $_2$	Implicit
24	OCLO	OClO	Implicit
25	CL2O2	Cl $_2$ O $_2$	Implicit
26	HCL	HCl	Implicit	wet
27	HOCL	HOCl	Implicit	wet
28	ClONO2	ClONO $_2$	Implicit	wet
29	BR	Br	Implicit
30	BRO	BrO	Implicit
31	HOBR	HOBr	Implicit	wet
32	HBR	HBr	Implicit	wet
33	BrONO 2	BrONO $_2$	Implicit	wet
34	BRCL	BrCl	Implicit

Table: Neutral Chemical Species (51 computed species + N $_2$ )

[mozart1a]

no.	Symbolic Name	Chemical Formula	Numerics	Deposition	Boundary Condition
35	CH4	CH $_4$	Implicit		vmr, airflux
36	CH3O2	CH $_3$ O $_2$	Implicit
37	CH3OOH	CH $_3$ OOH	Implicit	dry, wet
38	CH2O	CH $_2$ O	Implicit	dry, wet	flux
39	CO	CO	Explicit	dry	flux, ubvmr, airflux
40	CH3CL	CH $_3$ Cl	Explicit		vmr
41	CH3BR	CH $_3$ Br	Explicit		vmr
42	CFC11	CFCl $_3$	Explicit		vmr
43	CFC12	CF $_2$ Cl $_2$	Explicit		vmr
44	CFC113	CCl $_2$ FCClF $_2$	Explicit		vmr
45	HCFC22	CHClF $_2$	Explicit		vmr
46	CCL4	CCl $_4$	Explicit		vmr
47	CH3CCL3	CH $_3$ CCl $_3$	Explicit		vmr
48	CF2CLBR	CBr $_2$ F $_2$ (Halon-1211)	Explicit		vmr
49	CF3BR	CBrF $_3$ (Halon-1301)	Explicit		vmr
50	H2O	H $_2$ O	Explicit		flux
51	N2O	N $_2$ O	Explicit		vmr
52	CO2	CO $_2$	Explicit		vmr, ubvmr

Deposition:
wet = wet deposition included
dry = surface dry deposition included
If there is no designation in the deposition column, this species is not operated on by wet
or dry deposition algorthims.

Boundary Condition:
flux = flux lower boundary conditions
vmr = fixed volume mixing ratio (vmr) lower boundary condition
ubvmr = fixed vmr upper boundary condition
lflux = lightning emission included for this species
airflux= aircraft emissions included for this species
If there is no designation in the Boundary Conditions column, this species has a zero flux
boundary condition for the top and bottom of the model domain.

Table: (continued) Neutral Chemical Species (51 computed species + N $_2$ )

[mozart1b]

no.	Reactions	Comments
	Oxygen Reactions
1	O + O $_2$ + M $\rightarrow$ O $_3$ + M	JPL-06
2	O + O $_3$ $\rightarrow$ 2 O $_2$	JPL-06
3	O + O + M $\rightarrow$ O $_2$ + M	Smith and Robertson (2008)
4	O $_2$ ( $^1\Sigma$ ) + O $\rightarrow$ O $_2$ ( $^1\Delta$ ) + O	JPL-06
5	O $_2$ 1S + O $_2$ $\rightarrow$ O $_2$ ( $^1\Delta$ ) + O $_2$	JPL-06
6	O $_2$ ( $^1\Sigma$ ) + N $_2$ $\rightarrow$ O $_2$ ( $^1\Delta$ ) + N $_2$	JPL-06
7	O $_2$ ( $^1\Sigma$ ) + O $_3$ $\rightarrow$ O $_2$ ( $^1\Delta$ ) + O $_3$	JPL-06
8	O $_2$ ( $^1\Sigma$ ) + CO $_2$ $\rightarrow$ O $_2$ ( $^1\Delta$ ) + CO $_2$	JPL-06
9	O $_2$ ( $^1\Sigma$ ) $\rightarrow$ O $_2$	JPL-06
10	O $_2$ ( $^1\Delta$ ) + O $\rightarrow$ O $_2$ + O	JPL-06
11	O $_2$ ( $^1\Delta$ ) + O $_2$ $\rightarrow$ 2 O $_2$	JPL-06
12	O $_2$ ( $^1\Delta$ ) + N $_2$ $\rightarrow$ O $_2$ + N $_2$	JPL-06
13	O $_2$ ( $^1\Delta$ ) $\rightarrow$ O $_2$	JPL-06
14	O( $^1$ D) + N $_2$ $\rightarrow$ O + N $_2$	JPL-06
15	O( $^1$ D)+ O $_2$ $\rightarrow$ O + O $_2$ ( $^1\Sigma$ )	JPL-06
16	O( $^1$ D)+ O $_2$ $\rightarrow$ O + O $_2$	JPL-06
17	O( $^1$ D)+ H $_2$ O $\rightarrow$ 2 OH	JPL-06
18	O( $^1$ D) + N $_2$ O $\rightarrow$ 2 NO	JPL-06
19	O( $^1$ D) + N $_2$ O $\rightarrow$ N $_2$ + O $_2$	JPL-06
20	O( $^1$ D) + O $_3$ $\rightarrow$ 2 O $_2$	JPL-06
21	O( $^1$ D) + CFC11 $\rightarrow$ 3 Cl	JPL-06; Bloomfield (1994)
		for quenching of O( $^1$ D)
22	O( $^1$ D) + CFC12 $\rightarrow$ 2 Cl	JPL-06; Bloomfield (1994)
23	O( $^1$ D) + CFC113 $\rightarrow$ 3 Cl	JPL-06; Bloomfield (1994)
24	O( $^1$ D) + HCFC22 $\rightarrow$ Cl	JPL-06; Bloomfield (1994)
25	O( $^1$ D) + CCl $_4$ $\rightarrow$ 4 Cl	JPL-06
26	O( $^1$ D) + CH $_3$ Br $\rightarrow$ Br	JPL-06
27	O( $^1$ D) + CF $_2$ ClBr $\rightarrow$ Cl + Br	JPL-06
28	O( $^1$ D) + CF $_3$ Br $\rightarrow$ Br	JPL-06
29	O( $^1$ D) + CH $_4$ $\rightarrow$ CH $_3$ O $_2$ + OH	JPL-06
30	O( $^1$ D) + CH $_4$ $\rightarrow$ CH $_2$ O + H + HO $_2$	JPL-06
31	O( $^1$ D) + CH $_4$ $\rightarrow$ CH $_2$ O + H $_2$	JPL-06
32	O( $^1$ D) + H $_2$ $\rightarrow$ H + OH	JPL-06
33	O( $^1$ D) + HCl $\rightarrow$ Cl + OH	JPL-06
34	O( $^1$ D) + HBr $\rightarrow$ Br + OH	JPL-06

Table: Gas-phase Reactions.

no.	Reactions	Comments
	Nitrogen Radicals
35	N( $^2$ D) + O $_2$ $\rightarrow$ NO + O( $^1$ D)	$\qquad \qquad$ JPL-06 $\qquad \qquad$
36	N( $^2$ D) + O $\rightarrow$ N + O	JPL-06
37	N + O $_2$ $\rightarrow$ NO + O	JPL-06
38	N + NO $\rightarrow$ N $_2$ + O	JPL-06
39	N + NO $_2$ $\rightarrow$ N $_2$ O + O	JPL-06
40	NO + O + M $\rightarrow$ NO $_2$ + M	JPL-06
41	NO + HO $_2$ $\rightarrow$ NO $_2$ + OH	JPL-06
42	NO + O $_3$ $\rightarrow$ NO $_2$ + O $_2$	JPL-06
43	NO $_2$ + O $\rightarrow$ NO + O $_2$	JPL-06
44	NO $_2$ + O + M $\rightarrow$ NO $_3$ + M	JPL-06
45	NO $_2$ + O $_3$ $\rightarrow$ NO $_3$ + O $_2$	JPL-06
46	NO $_2$ + NO $_3$ + M $\rightarrow$ N $_2$ O5 + M	JPL-06
47	N $_2$ O $_5$ + M $\rightarrow$ NO $_2$ + NO $_3$ + M	JPL-06
48	NO $_2$ + OH + M $\rightarrow$ HNO $_3$ + M	JPL-06
49	HNO $_3$ + OH $\rightarrow$ NO $_3$ + H $_2$ O	JPL-06
50	NO $_2$ + HO $_2$ + M $\rightarrow$ HO $_2$ NO $_2$ + M	JPL-06
51	NO $_3$ + NO $\rightarrow$ 2 NO $_2$	JPL-06
52	NO $_3$ + O $\rightarrow$ NO $_2$ + O $_2$	$\qquad \qquad$ JPL-06 $\qquad \qquad$
53	NO $_3$ + OH $\rightarrow$ NO $_2$ + HO $_2$	JPL-06
54	NO $_3$ + HO $_2$ $\rightarrow$ NO $_2$ + OH + O $_2$	JPL-06
55	HO $_2$ NO $_2$ + OH $\rightarrow$ NO $_2$ + H $_2$ O + O $_2$	JPL-06
56	HO $_2$ NO $_2$ + M $\rightarrow$ HO $_2$ + NO $_2$ + M	JPL-06

Table: (continued) Gas-phase Reactions.

no.	Reactions	Comments
	Hydrogen Radicals
57	H + O $_2$ + M $\rightarrow$ HO $_2$ + M	$\qquad \qquad$ JPL-06 $\qquad \qquad$
58	H + O $_3$ + M $\rightarrow$ OH + O $_2$	JPL-06
59	H + HO $_2$ $\rightarrow$ 2 OH	JPL-06
60	H + HO $_2$ $\rightarrow$ H $_2$ + O $_2$	JPL-06
61	H + HO $_2$ $\rightarrow$ H $_2$ O + O	JPL-06
62	OH + O $\rightarrow$ H + O $_2$	JPL-06
63	OH + O $_3$ $\rightarrow$ HO $_2$ + O $_2$	JPL-06
64	OH + HO $_2$ $\rightarrow$ H $_2$ O + O $_2$	JPL-06
65	OH + OH $\rightarrow$ H $_2$ O + O	JPL-06
66	OH + OH + M $\rightarrow$ H $_2$ O $_2$ + M	JPL-06
67	OH + H $_2$ $\rightarrow$ H $_2$ O + H	JPL-06
68	OH + H $_2$ O $_2$ $\rightarrow$ H $_2$ O + HO $_2$	JPL-06
69	HO $_2$ + O $\rightarrow$ OH + O $_2$	JPL-06
70	HO $_2$ + O $_3$ $\rightarrow$ OH + 2O $_2$	JPL-06
71	HO $_2$ + HO $_2$ $\rightarrow$ H $_2$ O $_2$ + O $_2$	JPL-06
72	H $_2$ O $_2$ + O $\rightarrow$ OH + HO $_2$	JPL-06
	Chlorine Radicals
73	Cl + O $_3$ $\rightarrow$ ClO + O $_2$	JPL-06
74	Cl + H $_2$ $\rightarrow$ HCl + H	JPL-06
75	Cl + H $_2$ O $_2$ $\rightarrow$ HCl + HO $_2$	JPL-06
76	Cl + HO $_2$ $\rightarrow$ HCl + O $_2$	JPL-06
77	Cl + HO $_2$ $\rightarrow$ ClO + OH	JPL-06
78	Cl + CH $_2$ O $\rightarrow$ HCl + HO $_2$ + CO	JPL-06
79	Cl + CH $_4$ $\rightarrow$ CH $_3$ O $_2$ + HCl	JPL-06
80	ClO + O $\rightarrow$ Cl + O $_2$	JPL-06
81	ClO + OH $\rightarrow$ Cl + HO $_2$	JPL-06
82	ClO + OH $\rightarrow$ HCl + O $_2$	JPL-06
83	ClO + HO $_2$ $\rightarrow$ HOCl + O $_2$	JPL-06
84	ClO + NO $\rightarrow$ NO $_2$ + Cl	$\qquad \qquad$ JPL-06 $\qquad \qquad$
85	ClO + NO $_2$ + M $\rightarrow$ ClONO $_2$ + M	JPL-06

Table: (continued) Gas-phase Reactions.

no.	Reactions	Comments
	Chlorine Radicals Continued
86	ClO + ClO $\rightarrow$ 2 Cl + O $_2$	$\qquad \qquad$ JPL-06 $\qquad \qquad$
87	ClO + ClO $\rightarrow$ Cl2 + O $_2$	JPL-06
88	ClO + ClO $\rightarrow$ Cl + OClO	JPL-06
89	ClO + ClO + M $\rightarrow$ Cl $_2$ O $_2$ + M	JPL-06
90	Cl $_2$ O $_2$ + M $\rightarrow$ 2 ClO + M	JPL-06
91	HCl + OH $\rightarrow$ H $_2$ O + Cl	JPL-06
92	HCl + O $\rightarrow$ Cl + OH	JPL-06
93	HOCl + O $\rightarrow$ ClO + OH	JPL-06
94	HOCl + Cl $\rightarrow$ HCl + ClO	JPL-06
95	HOCl + OH $\rightarrow$ ClO + H $_2$ O	JPL-06
96	ClONO $_2$ + O $\rightarrow$ ClO + NO $_3$	JPL-06
97	ClONO $_2$ + OH $\rightarrow$ HOCl + NO $_3$	JPL-06
98	ClONO $_2$ + Cl $\rightarrow$ Cl $_2$ + NO $_3$	JPL-06
no.	Reactions	Comments
	Bromine Radicals
99	Br + O $_3$ $\rightarrow$ BrO + O $_2$	JPL-06
100	Br + HO $_2$ $\rightarrow$ HBr + O $_2$	JPL-06
101	Br + CH $_2$ O $\rightarrow$ HBr + HO $_2$ + CO	JPL-06
102	BrO + O $\rightarrow$ Br + O $_2$	JPL-06
103	BrO + OH $\rightarrow$ Br + HO $_2$	JPL-06
104	BrO + HO $_2$ $\rightarrow$ HOBr + O $_2$	JPL-06
105	BrO + NO $\rightarrow$ Br + NO $_2$	JPL-06
106	BrO + NO $_2$ + M $\rightarrow$ BrONO $_2$ + M	JPL-06
107	BrO + ClO $\rightarrow$ Br + OClO	JPL-06
108	BrO + ClO $\rightarrow$ Br + Cl + O $_2$	JPL-06
109	BrO + ClO $\rightarrow$ BrCl + O $_2$	JPL-06
110	BrO + BrO $\rightarrow$ 2 Br + O $_2$	JPL-06
111	HBr + OH $\rightarrow$ Br + H $_2$ O	JPL-06
112	HBr + O $\rightarrow$ Br + OH	JPL-06
113	HOBr + O $\rightarrow$ BrO + OH	JPL-06
114	BrONO $_2$ + O $\rightarrow$ BrO + NO $_3$	JPL-06

Table: (continued) Gas-phase Reactions.

no.	Reactions	Comments
	Halogen Radicals
115	CH $_3$ Cl + Cl $\rightarrow$ HO $_2$ + CO + 2HCl	JPL-06
116	CH $_3$ Cl + OH $\rightarrow$ Cl + H $_2$ O + HO $_2$	JPL-06
117	CH $_3$ CCl $_3$ + OH $\rightarrow$ 3 Cl + H $_2$ O	JPL-06
118	HCFC22 + OH $\rightarrow$ Cl + H $_2$ O + HO $_2$	JPL-06
119	CH $_3$ Br + OH $\rightarrow$ Br + H $_2$ O + HO $_2$	JPL-06
	CH $_4$ and Derivatives
120	CH $_4$ + OH $\rightarrow$ CH $_3$ O $_2$ + H $_2$ O	$\qquad \qquad$ JPL-06 $\qquad \qquad$
121	CH $_3$ O $_2$ + NO $\rightarrow$ CH $_2$ O + NO $_2$ + HO $_2$	JPL-06
122	CH $_3$ O $_2$ + HO $_2$ $\rightarrow$ CH $_3$ OOH + O $_2$	JPL-06
123	CH $_3$ OOH + OH $\rightarrow$ 0.7 CH $_3$ O $_2$ + 0.3 OH + 0.3 CH $_2$ O + H $_2$ O	JPL-06
124	CH $_2$ O + NO $_3$ $\rightarrow$ CO + HO $_2$ + HNO $_3$	JPL-06
125	CH $_2$ O + OH $\rightarrow$ CO + H $_2$ O + H	JPL-06
126	CH $_2$ O + O $\rightarrow$ OH + HO $_2$ + CO	JPL-06
127	CO + OH $\rightarrow$ H + CO $_2$	JPL-06

Table: (continued) Gas-phase Reactions.

no.	Reaction	Comments
	Sulfate Aerosol
1	N $_2$ O $_5$ + H $_2$ O $\rightarrow$ 2 HNO $_3$	JPL-06; f (sulfuric acid wt %)
2	ClONO $_2$ + H $_2$ O $\rightarrow$ HOCl + HNO $_3$	JPL-06; f (T, P, HCl, H $_2$ O, r)
3	BrONO $_2$ + H $_2$ O $\rightarrow$ HOBr + HNO $_3$	JPL-06; f (T, P, H $_2$ O, r)
4	ClONO $_2$ + HCl $\rightarrow$ Cl $_2$ + HNO $_3$	JPL-06; f (T, P, HCl, H $_2$ O, r)
5	HOCl + HCl $\rightarrow$ Cl $_2$ + H $_2$ O	JPL-06; f (T, P, HCl, HCl, H $_2$ O, r)
6	HOBr + HCl $\rightarrow$ BrCl + H $_2$ O	JPL-06; f (T, P, HCl, HOBr, H $_2$ O, r)
	NAT Aerosol
7	N $_2$ O $_5$ + H $_2$ O $\rightarrow$ 2 HNO $_3$	JPL-06; $\gamma = 4 \times 10^{-4}$
8	ClONO $_2$ + H $_2$ O $\rightarrow$ HOCl + HNO $_3$	JPL-06; $\gamma = 4 \times 10^{-3}$
9	ClONO $_2$ + HCl $\rightarrow$ Cl $_2$ + HNO $_3$	JPL-06; $\gamma = 0.2$
10	HCl + HCl $\rightarrow$ Cl $_2$ + H $_2$ O	JPL-06; $\gamma = 0.1$
11	BrONO $_2$ + H $_2$ O $\rightarrow$ HOBr + HNO $_3$	JPL-06; $\gamma = 0.3$
	Water-Ice Aerosol
12	N $_2$ O $_5$ + H $_2$ O $\rightarrow$ 2 HNO $_3$	JPL-06; $\gamma = 0.02$
13	ClONO $_2$ + H $_2$ O $\rightarrow$ HOCl + HNO $_3$	JPL-06; $\gamma = 0.3$
14	BrONO $_2$ + H $_2$ O $\rightarrow$ HOBr + HNO $_3$	JPL-06; $\gamma = 0.3$
15	ClONO $_2$ + HCl $\rightarrow$ Cl $_2$ + HNO $_3$	JPL-06; $\gamma = 0.3$
16	HOCl + HCl $\rightarrow$ Cl $_2$ + H $_2$ O	JPL-06; $\gamma = 0.2$
17	HOBr + HCl $\rightarrow$ BrCl + H $_2$ O	JPL-06; $\gamma = 0.3$

Table: Heterogeneous Reactions on liquid and solid aerosols.

no.	Reactants	Products	Comments
1	O $_2$ + h $\nu$	O + O( $^1$ D)	Ly- $\alpha$ : Chabrillat and Kockarts (1997, 1998)
			$\phi$ (Ly- $\alpha$ ): Lacoursiere et al. (1999)
			SRB: Koppers and Murtaugh (1996)
			For wavelength $\nu$ regions not Ly- $\alpha$ or SRB,
			$\sigma$ (120-205nm): Brasseur and Solomon (1986);
			$\sigma$ (205-240 nm): Yoshino et al. (1988)
2	O $_2$ + h $\nu$	2 O	see above
3	O $_3$ + h $\nu$	O( $^1$ D) + O $_2$	$\sigma$ (120-136.5nm): Tanaka et al. (1953);
			$\sigma$ (136.5-175nm): Ackerman (1971);
			$\sigma$ (175-847nm): WMO (1985); except for
			$\sigma$ (185-350nm): Molina and Molina (1986)
			$\phi$ ( $<$ 280nm): Marsh (1999)
			$\phi$ ( $>$ 280nm): JPL-06.
4	O $_3$ + h $\nu$	O + O $_2$	see above
5	N $_2$ O + h $\nu$	O( $^1$ D) + N $_2$	JPL-06
6	NO + h $\nu$	N + O	Minschwaner et al. (1993)
7	NO + h $\nu$	NO $^+$ + e
8	NO $_2$ + h $\nu$	NO + O	JPL-06
9	N $_2$ O $_5$ + h $\nu$	NO $_2$ + NO $_3$	JPL-06
10	N $_2$ O5 + h $\nu$	NO + O + NO $_3$	JPL-06
11	HNO $_3$ + h $\nu$	OH + NO $_2$	JPL-06
12	NO $_3$ + h $\nu$	NO $_2$ + O	JPL-06
13	NO $_3$ + h $\nu$	NO + O $_2$	JPL-06
14	HO $_2$ NO $_2$ + h $\nu$	OH + NO $_3$	JPL-06
15	HO $_2$ NO $_2$ + h $\nu$	NO $_2$ + HO $_2$	JPL-06
16	CH $_3$ OOH + h $\nu$	CH $_2$ O + H + OH	JPL-06
17	CH $_2$ O + h $\nu$	CO + 2 H	JPL-06
18	CH $_2$ O + h $\nu$	CO + H $_2$	JPL-06
19	H $_2$ O + h $\nu$	H + OH	$\phi$ (Ly- $\alpha$ ): Slanger et al. (1982);
			$\phi$ (105-145nm): Stief et al. (1975);
			$\phi$ ( $>$ 145): JPL-06
			$\phi$ (120-182nm): Yoshino et al. (1996);
			$\phi$ (183-194nm): Cantrell et al. (1997)

Table: Photolytic Reactions.

no.	Reactants	Products	Comments
20	H $_2$ O + h $\nu$	H $_2$ + O( $^1$ D)	(see above)
21	H $_2$ O + h $\nu$	H + 2 O	(see above)
22	H $_2$ O $_2$ + h $\nu$	2 OH	JPL-06
23	Cl $_2$ + h $\nu$	2 Cl	JPL-06
24	ClO + h $\nu$	Cl + O	JPL-06
25	OClO + h $\nu$	O + ClO	JPL-06
26	Cl $_2$ O $_2$ + h $\nu$	Cl + ClOO	Burkholder et al. (1990);
			Stimpfle et al. (2004)
27	HOCl + h $\nu$	Cl + OH	JPL-06
28	HCl + h $\nu$	Cl + H	JPL-06
29	ClONO $_2$ + h $\nu$	Cl + NO $_3$	JPL-06
30	ClONO $_2$ + h $\nu$	ClO + NO $_2$	JPL-06
31	BrCl + h $\nu$	Br + Cl	JPL-06
32	BrO + h $\nu$	Br + O	JPL-06
33	HOBr + h $\nu$	Br + OH	JPL-06
34	BrONO $_2$ + h $\nu$	Br + NO $_3$	JPL-06
35	BrONO $_2$ + h $\nu$	BrO + NO $_2$	JPL-06
36	CH $_3$ Cl + h $\nu$	Cl + CH $_3$ O $_2$	JPL-06
37	CCl $_4$ + h $\nu$	4 Cl	JPL-06
38	CH $_3$ CCl3 + h $\nu$	3 Cl	JPL-06
39	CFC11 + h $\nu$	3 Cl	JPL-06
40	CFC12 + h $\nu$	2 Cl	JPL-06
41	CFC113 + h $\nu$	3 Cl	JPL-06
42	HCFC22 + h $\nu$	Cl	JPL-06
43	CH $_3$ Br + h $\nu$	Br + CH $_3$ O $_2$	JPL-06
44	CF $_3$ Br + h $\nu$	Br	JPL-06
45	CF $_2$ ClBr + h $\nu$	Br + Cl	JPL-06
46	CO $_2$ + h $\nu$	CO + O	$\sigma$ (120-167): Nakata, et al. (1965);
			$\sigma$ (167-199): Huffman (1971)
47	CH $_4$ + h $\nu$	H + CH $_3$ O $_2$	$\sigma$ : JPL-06;
			based on Brownsword et al. (1997)
48	CH $_4$ + h $\nu$	H $_2$ + 0.18 CH $_2$ O + 0.18 O
		0.44 CO $_2$ + 0.44 H $_2$	see above
		0.38 CO + 0.05 H $_2$ O

Table: (continued) Photolytic Reactions.

Reaction	$\Delta H$ (kJ mol $^{-1}$ )
O $^+$ + O $_2$ $\rightarrow$ O $_2^+$ + O	150.11
O $^+$ + N $_2$ $\rightarrow$ NO $^+$ + N	105.04
N $_2^+$ + O $\rightarrow$ NO $^+$ + N( $^2$ D)	67.53
O $_2^+$ + N $\rightarrow$ NO $^+$ + O	406.16
O $_2^+$ + NO $\rightarrow$ NO $^+$ + O $_2$	271.38
N $^+$ + O $_2$ $\rightarrow$ O $_2^+$ + N	239.84
N $^+$ + O $_2$ $\rightarrow$ NO $^+$ + O	646.28
N $^+$ + O $\rightarrow$ O $^+$ + N	95.55
N $_2^+$ + O $_2$ $\rightarrow$ O $_2^+$ + N $_2$	339.59
O $_2^+$ + N $_2$ $\rightarrow$ NO $^+$ + NO	–
N $_2^+$ + O $\rightarrow$ O $^+$ + N $_2$	–
NO $^+$ + e $\rightarrow$ 0.2N + 0.8N( $^2$ D) + O	82.389
O $_2^+$ + e $\rightarrow$ 1.15O + 0.85O( $^1$ D)	508.95
N $_2^+$ + e $\rightarrow$ 1.1N + 0.9N( $^2$ D)	354.83

Table: Ion-neutral and recombination reactions and exothermicities.

[ionrxntab]

O + h $\nu$ $\rightarrow$ O $^+$ + e

O + e $^*$ $\rightarrow$ O $^+$ + e + e $^*$

N + hv $\rightarrow$ N $^+$ + e

O $_2$ + h $\nu$ $\rightarrow$ O $_2^+$ + e

O $_2$ + e $^*$ $\rightarrow$ O $_2^+$ + e + e $^*$

O $_2$ + h $\nu$ $\rightarrow$ O + O $^+$ + e

O $_2$ + e $^*$ $\rightarrow$ O + O $^+$ + e + e $^*$

N $_2$ + h $\nu$ $\rightarrow$ N $_2^+$ + e

N $_2$ + e $^*$ $\rightarrow$ N $_2^+$ + e + e $^*$

N $_2$ + h $\nu$ $\rightarrow$ N + N $^+$ + e

N $_2$ + e $^*$ $\rightarrow$ N + N $^+$ + e + e $^*$

N $_2$ + h $\nu$ $\rightarrow$ N( $^2$ D) + N $^+$ + e

N $_2$ + e $^*$ $\rightarrow$ N( $^2$ D) + N $^+$ + e + e $^*$

Table: Ionization reactions.

wavelength interval	$F_i^0$	$R_i$
nm	ph cm $^{-2}$ s $^{-1}$
0.05 - 0.4	5.010e+01	6.240e-01
0.4 - 0.8	1.000e+04	3.710e-01
0.8 - 1.8	2.000e+06	2.000e-01
1.8 - 3.2	2.850e+07	6.247e-02
3.2 - 7.0	5.326e+08	1.343e-02
7.0 - 15.5	1.270e+09	9.182e-03
15.5 - 22.4	5.612e+09	1.433e-02
22.4 - 29.0	4.342e+09	2.575e-02
29.0 - 32.0	8.380e+09	7.059e-03
32.0 - 54.0	2.861e+09	1.458e-02
54.0 - 65.0	4.830e+09	5.857e-03
65.0 - 79.8	1.459e+09	5.719e-03
65.0 - 79.8	1.142e+09	3.680e-03
79.8 - 91.3	2.364e+09	5.310e-03
79.8 - 91.3	3.655e+09	5.261e-03
79.8 - 91.3	8.448e+08	5.437e-03
91.3 - 97.5	3.818e+08	4.915e-03
91.3 - 97.5	1.028e+09	4.955e-03
91.3 - 97.5	7.156e+08	4.422e-03
97.5 - 98.7	4.482e+09	3.950e-03
98.7 - 102.7	4.419e+09	5.021e-03
102.7 - 105.0	4.235e+09	4.825e-03
105.0 - 121.0	2.273e+10	3.383e-03

Table: EUVAC model parameters.

[euvactab]

7.6.8. Electric Field¶

The global electric field is based on a composite of two empirical models for the different latitude regions: at high latitude the Weimer95 model (???), and at low- and midlatitude the Scherliess model (???). In the following the different models are described since the model is not published to date.

7.6.8.1. Low- and midlatitude electric potential model¶

The low- and mid latitude electric field model was developed by Ldger Scherliess (???). It’s based on Incoherent Scatter Radar data (ISR) from Jicamarca, Arecibo, Saint Santin, Millstone Hill, and the MU radar in Shigaraki. The electric field is calculated for a given year, season, UT, $S_a$ , local time, and with longitudinal/latitudinal variation. The empirical model is constructed from a model for low solar flux ( $S_a = 90$ ) and a high solar flux model ( $S_a = 180$ ). The global electric potential is expressed according to (???) by

$\begin{split} \Phi(d,T,t,\lambda) = \sum_{k=0}^2 \sum_{l=-2}^2 \sum_{m=-n}^n \sum_{n=1}^{12} & A_{klmn} P_n^m(sin \lambda) f_m(\frac{2 \Pi t}{24}) \\ & f_l(\frac{2 \Pi T}{24}) f_{-k}(\frac{2 \Pi (d + 9)}{365.24}) \end{split}$

with

$\begin{aligned} f_m(\phi) & = \sqrt{2} \sin(m \phi) \quad & m > 0 \\ f_m(\phi) & = 1 \quad & m = 0 \\ f_m(\phi) & = \sqrt{2} \cos(m \phi) \quad & m < 0\end{aligned}$

the day of the year is denoted by $d$ , universal time by $T$ , magnetic local time by $t$ , and geomagnetic latitude $\lambda$ . The values of $d$ , $T$ , and $t$ are expressed as angles between 0 and 2 $\Pi$ . $P_n^m$ are fully normalized Legendre polynomials. Due to the assumption that the geomagnetic field lines are highly conducting, the $n+m$ odd coefficients are set to zero to get a symmetrical electric potential about the magnetic equator. The coefficients $A_{klmn}$ are found by a least–square fit for low and high solar flux. The solar cycle dependence is introduced by inter- and extrapolation of the sets of coefficients $A_{klmn}^{low}$ for $S_a = 90$ and $A_{klmn}^{high}$ for $S_a = 180$ .

$\begin{aligned} A_{klmn} =A^{low}_{klmn} + S_{aM}[A^{high}_{klmn} - A^{low}_{klmn}] \end{aligned}$

with

$\begin{aligned} S_{aM} & = \frac{arctan[(S_a-65)^2/ 90^2]- a_{90}}{a_{180}-a_{90}} \\ a_{90} & = arctan[(90-65)^2/90^2]\\ a_{180} & = arctan[(180-65)^2/90^2]\end{aligned}$

We are using the daily $F_{10.7}$ number for $S_a$ . $S_{aM}$ levels off at high and low solar flux numbers, and therefore the model does not predict unrealistic high or low electric potential values.

The geomagnetic field is described by modified apex coordinates (???) which already take into account the distortion of the magnetic field. Modified apex coordinates have a reference height associated with them, which in our case is set to 130 km. The electric field $\mathbf{E}$ and the electromagnetic drift velocity $\mathbf{v}_E$ can be expressed by quantities mapped to the reference height, e.g. by $E_{d1}$ , $E_{d2}$ and $v_{e1}$ , $v_{e2}$ . These quantities are not actual electric field or electromagnetic drift velocity components, but rather the representation of the electric field or electromagnetic drift velocities by being constant along the geomagnetic field line. The fields in an arbitrary direction $\mathbf{I}$ can be expressed by

$\begin{aligned} \mathbf{I} \cdot \mathbf{E} &=\mathbf{I} \cdot \mathbf{d}_1 E_{d1} + \mathbf{I} \cdot \mathbf{d}_2 E_{d2} \\ \mathbf{I} \cdot \mathbf{v}_E &= \mathbf{I} \cdot \mathbf{e}_1 v_{e1} + \mathbf{I} \cdot \mathbf{e}_2 v_{e2}\end{aligned}$

The basis vector $\mathbf{d}_1$ and $\mathbf{e}_1$ are in more–or–less magnetic eastward direction and $\mathbf{d}_2$ and $\mathbf{e}_2$ in downward/ equatorward direction. The base vectors vary with height, $\mathbf{d}_i$ is decreasing and $\mathbf{e}_i$ increasing with altitude. Therefore when the base vectors are applied to the mapped field at the reference height, e.g. $E_{d1}$ , $E_{d2}$ and $v_{e1}$ , $v_{e2}$ , they already take into account the height and directional variation of the corresponding quantity. Note that the modified apex coordinates are using the International Geomagnetic Reference Field (IGRF), and in the WACCM4 code the IGRF is only defined between the years 1900 and 2000. The description of the IGRF can be updated every 5 years to be extended in time.

7.6.8.2. High–latitude electric potential model¶

The high–latitude electric potential model from Weimer (???) is used. The model is based on spherical harmonic coefficients that were derived by least square fitting of measurements from the Dynamics Explorer 2 (DE2) satellite. The variation of the spherical harmonic coefficients with the interplanetary magnetic field (IMF) clock angle, IMF strength, solar wind velocity and season can be reproduced by a combination of Fourier series and multiple linear regression formula. The final model varies with magnetic latitude, magnetic local time, season, IMF strength and direction, and solar wind velocity. For our purpose we have set the solar wind speed to a constant value of 400 km/s and only consider the effects of IMF $B_z$ ( $B_y=0$ ). Since the IMF conditions are not known all the time, we developed an empirical relation between $B_z$ and the $K_p$ index and the solar flux number $S_a$ . Both, the $K_p$ index and the daily solar flux number $F_{10.7}$ , are known in the WACCM4 model.

$\begin{split} B_z (K_p, F_{10.7}) = & - 0.085 K_p^2 - 0.08104 K_p + 0.4337 + \\ & 0.00794 F_{10.7} - 0.00219 K_p F_{10.7} \end{split}$

Note that the Weimer model uses an average year of 365.24 days/year and an average month of 30.6001 days/month. The boundary of the Weimer model is at $46^o$ magnetic latitude. The model was developed for an averaged northern and southern hemisphere. The $B_y$ value and the season are reversed to get the values for the other hemisphere.

7.6.8.3. Combing low–/ mid–latitude with the high latitude electric potential¶

After the low/mid–latitude electric potential $\Phi_{mid}$ and the high latitude potential $\Phi_{hgh}$ are calculated, both patterns are combined to be smooth at the boundary. The boundary between high and mid latitude $\lambda_{bnd}$ is defined to lie where the electric field magnitude $E$ from $\Phi_{hgh}$ equals 15 mV/m. After finding the longitudinal variation of the high latitude boundary $\lambda_{bnd}$ , it’s shifted halfway towards $54^o$ magnetic latitude. The width of the transition zone $2 \Delta \lambda_{trs}$ from high to mid latitude varies with magnetic local time. First, the high and mid latitude electric potential are adjusted by a constant factor such that the average for the high and mid latitude electric potential along the boundary $\lambda_{bnd}$ are the same. The combined electric potential $\Phi$ is defined by

$\Phi = \begin{cases} \Phi_{mid} \quad & | \lambda | < \lambda_{bnd}-\Delta \lambda_{trs} \\ \Phi_{hgh} \quad & | \lambda | > \lambda_{bnd}+\Delta \lambda_{trs} \\ F_{int}(\Phi_{mid},\Phi_{hgh}) \quad & \lambda_{bnd}-\Delta \lambda_{trs} \leq | \lambda | \leq \lambda_{bnd}+\Delta \lambda_{trs} \end{cases}$

with

$\begin{split} F_{int}(\Phi_{mid},\Phi_{hgh}) = \frac{1}{3} \frac{1}{2 \Delta \lambda_{trs}}[ & \left\{ \Phi_{mid}(\phi,\lambda_{bnd}-\Delta \lambda_{trs}) + 2 \Phi_{mid}(\phi,\lambda) \right\} \\ & \left\{\lambda_{bnd}-|\lambda|+\Delta \lambda_{trs} \right\} + ( \Phi_{hgh}(\phi,\lambda_{bnd}+\Delta \lambda_{trs}) + \\ & 2 \Phi_{hgh}(\phi,\lambda) ) \left\{-\lambda_{bnd}+|\lambda|+\Delta \lambda_{trs}\right\} ] \end{split}$

7.6.8.4. Calculation of electric field¶

The electric field can be derived from the electric potential by

$\begin{aligned} \mathbf{E} = - \nabla \Phi\end{aligned}$

The more-or-less magnetic eastward electric field component $E_{d1}$ and the in general downward/ equatorward $E_{d2}$ component are calculated. These components are constant along the magnetic field line. They are calculated at a reference height $h_r= \; 130$ km with $R = R_{earth}+ h_r$ . The electric field does not vary much with altitude, and therefore we assume in the code that the electric field is constant in height.

$\begin{aligned} E_{d1} &= -\frac{1}{R cos \lambda}\frac{\partial \Phi}{\partial \phi} \\ E_{d2} &= \frac{1}{R \sin I}\frac{\partial \Phi}{\partial \lambda}\end{aligned}$

with $\sin I = 2\sin \lambda [4-3\cos^2\lambda]^{0.5}$ .

7.6.8.5. Calculation of electrodynamic drift velocity¶

The electric field is calculated on a $2^o \; \times \; 2^o$ degree geomagnetic grid with the magnetic longitude represented by the magnetic local time (MLT) from 0 MLT to 24 MLT. Therefore, the magnetic local time of the geographic longitudes of the WACCM4 grid has to be determined first to map from the geomagnetic to the geographic WACCM4 grid. The magnetic local time is calculated by using the location of the geomagnetic dipole North pole, the location of the subsolar point, and the apex longitude of the geographic WACCM4 grid point. A bilinear interpolation is used for the mapping. Note that every processor calculates the global electric field, which is computationally inexpensive. Otherwise, to calculate the electric field some communication between the different processors would be necessary to get the spatial derivatives. The mapped electric field is rotated into the geographic direction by

$\begin{aligned} \mathbf{E} &= \mathbf{d}_1 E_{d1} + \mathbf{d}_2 E_{d2}\end{aligned}$

with the components of $\mathbf{E}$ being the geographic eastward, westward and upward electric field. At high altitudes the ion–neutral collision frequency $\nu_{in}$ is small in relation to the angular gyrofrequency of the ions $\Omega_i$ ( $\nu_{in} \ll \Omega_i$ ), and the electron–neutral collision frequency $\nu_{en}$ is much smaller than the angular gyrofrequency of the electrons $\Omega_e$ ( $\nu_{en} \ll \Omega_e$ ), due to the decrease in neutral density with increasing altitude. Therefore, the ion drift $\mathbf{v}_{i\bot}$ perpendicular to the geomagnetic field can be simplified by the electrodynamic drift velocity $\mathbf{v}_E$

$\begin{aligned} \mathbf{v}_{i\bot} \approx \mathbf{v}_E = \frac{\mathbf{E} \times \mathbf{B}_o}{\mathbf{B}_o^2}\end{aligned}$

with $\mathbf{B}_o$ the geomagnetic main field from IGRF.

7.6.8.6. Ion drag calculation¶

The following is written according to the source code. Two subroutines iondrag_calc exist in the code, one uses the calculated ion drag coefficients if WACCM_MOZART is used, and the other one uses look-up tables for the ion drag coefficients $\lambda_1$ and $\lambda_2$ . It is assumed that the electron $T_e$ and ion $T_i$ temperature is equal to the neutral temperature $T_n$ .

$\begin{aligned} T_i = T_e = T_n\end{aligned}$

The dip angle I of the geomagnetic field is calculated by

$\begin{aligned} I = \arctan \frac{B_z}{\sqrt{B_{north}^2+B_{east}^2}}\end{aligned}$

with a minimum dip angle $|I| \geq 0.17$ . The declination is

$\begin{aligned} D = \arctan \frac{B_{east}}{B_{north}}\end{aligned}$

The magnetic field component $B_z, B_{east}, B_{north}$ are determined from the International Geomagnetic Reference Field (IGRF). The collision frequencies $\nu$ in units of $s^{-1}$ are determined by, e.g. (???)

$\begin{aligned} \frac{1}{N_{O_2}} \nu_{O_2^+ - O_2} &= 2.59\times 10^{-11}\sqrt{\frac{T_i + T_e}{2}}\left[ 1-0.73 log_{10}\sqrt{\frac{T_i + T_e}{2}}\right]^2 \\ \frac{1}{N_{O_2}} \nu_{O^+ - O_2} &=6.64\times 10^{-10} \\ \frac{1}{N_{O_2}} \nu_{NO^+ - O_2} &=4.27\times 10^{-10} \\ \frac{1}{N_{O}} \nu_{O^+ - O} &=3.67\times 10^{-11}\sqrt{\frac{T_i + T_e}{2}}\left[ 1-0.064 log_{10}\sqrt{\frac{T_i + T_e}{2}}\right]^2 f_{cor} \\ \frac{1}{N_{O}} \nu_{NO^+ - O} &=2.44\times 10^{-10} \\ \frac{1}{N_{O}} \nu_{O_2^+ - O} &=2.31\times 10^{-10} \\ \frac{1}{N_{N_2}} \nu_{O_2^+ - N_2}&=4.13\times 10^{-10} \\ \frac{1}{N_{N_2}} \nu_{NO^+ - N_2} &=4.34\times 10^{-10} \\ \frac{1}{N_{N_2}} \nu_{O^+ - N_2} &=6.82\times 10^{-10}\end{aligned}$

with $N_n$ the number density for the neutral n in units of $1/cm^3$ , and the temperature in Kelvins. The collisions frequencies for $\nu_{O_2^+ - O_2}$ and $\nu_{O^+ - O}$ are resonant, all other are nonresonant. The arbitrary correction factor $f_{cor}$ multiplies the $\nu_{O^+ - O}$ collision frequency and is set to $f_{cor} =1.5$ which has been found to improve agreement between calculated and observed winds and electron densities in the upper thermosphere in other models. The mean mass $\overline{m}_{mid}$ [g/mole] at the midpoints of the height level is calculated in the Mozart module. The number densities $[1/cm^3]$ are

$\begin{aligned} N_{O_2} &= \frac{N \overline{m}_{mid} mmr_{O_2}}{m_{O_2}} \\ N_{O} &= \frac{N \overline{m}_{mid} mmr_{O}}{m_{O}} \\ N_{N_2} &= \frac{N \overline{m}_{mid} mmr_{N_2}}{m_{N_2}} \\ N_{O_2^+}&= \frac{N \overline{m}_{mid} mmr_{O_2^+}}{m_{O_2^+}} \\ N_{O^+} &= \frac{N \overline{m}_{mid} mmr_{O^+}}{m_{O^+}} \\ N_{e} &= \frac{N \overline{m}_{mid} mmr_e }{m_{e}}\end{aligned}$

with $mmr$ the mass mixing ratio, and $N$ the total number density in units of $1/cm^3$ . The pressure $[dyne/cm^2]$ and the mean mass at the midpoint $\overline{m}_{mid}$ in units of $g/mole$ are

$\begin{aligned} p = 10 \; p_{mid} \\ N \overline{m}_{mid} = \frac{p \; \overline{m}}{k_B T_n}\end{aligned}$

with the factor 10 to convert from [Pa] to $[dyne/cm^2]$ , and $k_B$ the Boltzmann constant. The collision frequencies are

$\begin{aligned} \nu_{O_2^+} &= \nu_{O_2^+ - O_2} + \nu_{O_2^+ - O} + \nu_{O_2^+ - N_2} \\ \nu_{O^+} &= \nu_{O^+ - O_2} + \nu_{O^+ - O} + \nu_{O^+ - N_2} \\ \nu_{NO^+} &= \nu_{NO^+ - O_2} + \nu_{NO^+ - O} + \nu_{NO^+ - N_2} \\ \begin{split} \nu_{en} &= 2.33\times10^{-11} N_{N_2} T_e (1-1.21 \times 10^{-4} T_e) + \\ & 1.82 \times 10^{-10} N_{O_2} \sqrt{T_e} (1 + 3.6 \times 10^{-2} \sqrt{T_e}) + \\ & 8.9 \times 10^{-11} N_O \sqrt{T_e} (1 + 5.7 \times 10^{-4} T_e) \end{split}\end{aligned}$

The ratios $r$ between collision frequency $\nu$ and gyro frequency $\Omega$ are

$\begin{aligned} r_{O_2^+} &= \frac{\nu_{O_2^+}}{\Omega_{O_2^+}}\\ r_{O^+} &= \frac{\nu_{O^+}}{\Omega_{O^+}}\\ r_{NO^+} &= \frac{\nu_{NO^+}}{\Omega_{NO^+}}\\ r_{e} &= \frac{\nu_{en}}{\Omega_{e}}\end{aligned}$

with the gyro frequency for ions $\Omega_i = e B/m_i$ and for electrons $\Omega_e=eB/m_e$ . The Pedersen conductivity [ S/m] is

$\begin{split} \sigma_P = &\frac{e}{B} [ N_{O^+} \frac{r_{O^+}} {1+r_{O^+}^2 } + N_{O_2^+} \frac{r_{O_2^+}} {1+r_{O_2^+}^2 } + \\ & N_{NO^+} \frac{r_{NO^+}} {1+r_{NO^+}^2 } + N_{e} \frac{r_e} {1+r_e^2 } ] \end{split}$

The Hall conductivity [S/m] is

$\begin{split} \sigma_H = &\frac{e}{B} [ -N_{O^+} \frac{1} {1+r_{O^+}^2 } - N_{O_2^+} \frac{1} {1+r_{O_2^+}^2 } - \\ & N_{NO^+} \frac{1} {1+r_{NO^+}^2 }+ N_{e} \frac{1} {1+r_{e}^2 } ] \end{split}$

The ion drag coefficients are

$\begin{aligned} \lambda_1 &= \frac{\sigma_P B^2}{\rho} \\ \lambda_2 &= \frac{\sigma_H B^2}{\rho}\end{aligned}$

with $\rho= N \frac{\overline{m}}{N_A}$ , and $N_A$ the Avagadro number. The ion drag tensor in magnetic direction $\underline{\lambda}^{mag}$ is

$\begin{gathered} \underline{\lambda}^{mag}= \begin{pmatrix} \lambda_{xx}^{mag} & \lambda_{xy}^{mag} \\ \lambda_{yx}^{mag} & \lambda_{yy}^{mag} \end{pmatrix} = \begin{pmatrix} \lambda_1 & \lambda_{2}sin I\\ -\lambda_2 sin I & \lambda_{1} sin^2 I \end{pmatrix}\end{gathered}$

with the x–direction in magnetic east, and y–direction magnetic north in the both hemispheres. The ion drag tensor can be rotated in geographic direction by using the rotation matrix $\mathbf{R}$

$\begin{gathered} \mathbf{R} = \begin{pmatrix} \cos D & \sin D\\ -\sin D & \cos D \end{pmatrix}\end{gathered}$

Applying the rotation to the ion drag tensor $\mathbf{R}\underline{\lambda}^{mag}\mathbf{R}^{-1}$ leads to

$\begin{gathered} \Lambda = \begin{pmatrix} \lambda_{xx} & \lambda_{xy} \\ \lambda_{yx} & \lambda_{yy} \end{pmatrix} = \\ \begin{pmatrix} \lambda_{xx}^{mag} cos^2 D + \lambda_{yy}^{mag} sin^2 D & \lambda_{xy}^{mag} + (\lambda_{yy}^{mag}- \lambda_{xx}^{mag}) \sin D \cos D \\ \lambda_{yx}^{mag} + (\lambda_{yy}^{mag}- \lambda_{xx}^{mag}) \sin D \cos D & \lambda_{yy}^{mag} \cos^2 D + \lambda_{xx}^{mag} \sin^2 D \end{pmatrix}\end{gathered}$

The ion drag acceleration $\mathbf{a}_i$ due to the Ampère force is

$\begin{aligned} \mathbf{a}_i = \frac{\mathbf{J}\times \mathbf{B}}{\rho} = \lambda_1 (\mathbf{v}_E - \mathbf{u}_{n\perp}) + \lambda_2 \mathbf{\hat{b}}\times (\mathbf{v}_E - \mathbf{u}_{n\perp})\end{aligned}$

with $\mathbf{u}_{n\perp}$ the neutral wind velocity perpendicular to the geomagnetic field and $\mathbf{\hat{b}}$ the unit vector of the geomagnetic field. The tendencies on the neutral wind are calculated by

$\begin{aligned} \frac{\partial \mathbf{v}_{En}}{\partial t} = -\Lambda \mathbf{v}_{En}\end{aligned}$

For stability an implicit scheme is used with

$\begin{aligned} \frac{\mathbf{v}_{En}(t+\Delta t) - \mathbf{v}_{En}(t)}{\Delta t} = -\Lambda \mathbf{v}_{En}(t+\Delta t)\end{aligned}$

which leads to

$\begin{aligned} (\frac{1}{\Delta t} I + \Lambda)\mathbf{v}_{En}(t+\Delta t) = \frac{1}{\Delta t}\mathbf{v}_{En}(t)\end{aligned}$

with $I$ the unit matrix. Solving for $\mathbf{v}_{En}(t+\Delta t)$ gives

$\begin{aligned} \mathbf{v}_{En}(t+\Delta t) = \frac{1}{\Delta t}(\frac{1}{\Delta t} I + \Lambda)^{-1}\mathbf{v}_{En}(t)\end{aligned}$

The tendencies are determined by

$\begin{aligned} \frac{\partial \mathbf{v}_{En}}{\partial t} = \frac{\mathbf{v}_{En}(t+\Delta t) - \mathbf{v}_{En}(t)}{\Delta t} = \frac{1}{\Delta t}[ \frac{1}{\Delta t}(\frac{1}{\Delta t} I + \Lambda)^{-1}-1] \mathbf{v}_{En}(t)\end{aligned}$

The tensor $\frac{1}{\Delta t} I + \Lambda$ is

$\begin{gathered} \begin{pmatrix} \lambda_{11}^* & \lambda_{12}^* \\ \lambda_{21}^* & \lambda_{22}^* \end{pmatrix} = \begin{pmatrix} \frac{1}{\Delta t} + \lambda_{xx}& \lambda_{xy} \\ \lambda_{yx} & \frac{1}{\Delta t } + \lambda_{yy} \end{pmatrix}\end{gathered}$

$\begin{aligned} \frac{Det}{\Delta t} = \frac{1}{\Delta t} \frac{1}{\lambda_{11}^* \lambda_{22}^* - \lambda_{12}^* \lambda_{21}^*}\end{aligned}$

The tendencies applied to the neutral winds with $\mathbf{v}_{En}=(u_E- u_n, v_E - v_n)$ gives

$\begin{aligned} d_t u_i =& \frac{1}{\Delta t} \left[\frac{Det}{\Delta t} \left( \lambda_{12}^* (v_E - v_n) - \lambda_{22}^* (u_E- u_n) \right) + u_E - u_n \right] \\ d_t v_i =& \frac{1}{\Delta t} \left[ \frac{Det}{\Delta t} \left( \lambda_{21}^* (u_E -u_n) - \lambda_{11}^* (v_E- v_n) \right) + v_E - v_n \right]\end{aligned}$

The electromagnetic energy transfer to the ionosphere is

$\begin{aligned} \mathbf{J}\cdot\mathbf{E} = \mathbf{J}\cdot\mathbf{E}' + \mathbf{u}_n \cdot \mathbf{J}\times\mathbf{B}\end{aligned}$

The first term on the right hand side denotes the Joule heating, which is the electromagnetic energy transfer rate in the frame of reference of the neutral wind. The second term represents the generation of kinetic energy due to the Ampère force. Since the electric field is small along the magnetic field line, we consider only the perpendicular component to the magnetic field of the Joule heating $\mathbf{J}_{\perp}\cdot\mathbf{E}'$ . The electric field in the frame of the neutral wind $\mathbf{u}$ can be written as

$\begin{aligned} \mathbf{E}' = \mathbf{E} + \mathbf{u}\times \mathbf{B}\end{aligned}$

The Joule heating can be expressed by

$\begin{aligned} \mathbf{J}_{\perp}\cdot\mathbf{E}' = \sigma_P \mathbf{E}'^2\end{aligned}$

with

$\begin{aligned} \mathbf{E}'^2 = B^2 (\frac{\mathbf{E}\times\mathbf{B}}{B^2} - \mathbf{u}_{\perp})^2\end{aligned}$

and $\frac{\mathbf{E}\times\mathbf{B}}{B^2}$ the electromagnetic drift velocity $\mathbf{v}_E$ with the components $u_E$ and $v_E$ . The Joule heating $Q_J$ is

$\begin{aligned} Q_J = (u_E - u_n)^2 \lambda_{xx} + (u_E - u_n)(v_E - v_n) (\lambda_{xy} - \lambda_{yx}) _+ (v_E - v_n)^2 \lambda_{yy}\end{aligned}$

Note, that the vertical velocity components are not taken into account here.

7.6.9. Boundary Conditions¶

The upper boundary conditions for momentum and for most constituents are the usual zero flux conditions used in CAM4. However, in the energy budget of the thermosphere, much of the SW radiation at wavelengths $<$ 120 nm is absorbed above 145 km (the upper boundary of the model), where LW radiation is very inefficient. This energy is transported downward by molecular diffusion to below 120 km, where it can be dissipated more efficiently by LW emission. Imposing a zero flux upper boundary condition on heat omits a major term in the heat budget and causes the lower thermosphere to be much too cold. Instead, we use the Mass Spectrometer-Incoherent Scatter (MSIS) model (???; ???) to specify the temperature at the top boundary as a function of season and phase of the solar cycle. The version of the MSIS model used in is NRLMSISE-00 [see http://uap-www.nrl.navy.mil/models_web/msis/msis_home.htm].

For chemical constituents, surface mixing ratios of CH $_4$ , N $_2$ O, CO $_2$ , H $_2$ , CFC-11, CFC-12, CFC-113, HCFC-22, H-1211, H-1301, CCl $_4$ , CH $_3$ CCH $_3$ , CH $_3$ Cl, and CH $_3$ Br are specified from observations. The model accounts for surface emissions of NO $_{\rm X}$ and CO based on the emission inventories described in (???). The NO $_{\rm X}$ source from lightning is distributed according to the location of convective clouds based on (???) and (???), with a vertical profile following (???). Aircraft emissions of NO $_{\rm X}$ and CO are included in the model and based on (???).

At the upper boundary, a zero-flux upper boundary condition is used for most species whose mixing ratio is negligible in the lower thermosphere, while mixing ratios of other species are specified from a variety of sources. The MSIS model is used to specify the mixing ratios of O, O $_2$ , H, and N; as in the case of temperature, the MSIS model returns values of these constituents as functions of season and phase of the solar cycle. CO and CO $_2$ are specified at the upper boundary using output from the TIME-GCM (???). NO is specified using data from the Student Nitric Oxide Explorer (SNOE) satellite (???), which has been parameterized as a function of latitude, season, and phase of the solar cycle in the Nitric Oxide Empirical Model (NOEM) of (???). Finally, a global-mean value (typical of the sunlit lower thermosphere) is specified for species such as H $_2$ O, whose abundance near the top of the model is very small under sunlit conditions, but which can be rapidly transported upward by diffusive separation in polar night (since they are lighter than the background atmosphere). In these cases, a zero-flux boundary condition leads to unrealistically large mixing ratios at the model top in polar night.

7. Extensions to CAM¶

7.1. Introduction¶

7.2. Chemistry¶

7.2.1. Chemistry Schemes¶

7.2.1.1. The Bulk Aerosol Model¶

7.2.1.2. CAM-Chem using the Modal Aerosol Model¶

7.2.1.3. Trop MOZART Chemistry¶

7.2.1.4. Trop-Strat MOZART Chemistry¶

7.2.1.5. SOA calculation in CAM-Chem¶

7.2.2. Emissions¶

7.2.2.1. Emissions of Trace Gases¶

7.2.2.2. Biogenic emissions¶

7.2.3. Boundary conditions¶

7.2.3.1. Lower boundary conditions¶

7.2.3.2. Specified stratospheric distributions¶

7.2.3.3. Upper boundary condition¶

7.2.4. Lightning¶

7.2.5. Dry deposition¶

7.2.6. Wet removal¶

7.3. Photolytic Approach (Neutral Species)¶

7.4. Numerical Solution Approach¶

7.5. Superfast Chemistry¶

7.5.1. Chemical mechanism¶

7.5.2. Emissions for CAM4 superfast chemistry¶

7.5.3. LINOZ¶

7.5.4. Parameterized PSC ozone loss¶

7.6. Physical Parameterizations¶

7.6.1. Domain and Resolution¶

7.6.2. Molecular Diffusion and Constituent Separation¶

7.6.2.1. Molecular viscosity and diffusivity¶

7.6.2.2. Diffusive separation velocities¶

7.6.2.3. Discretization of the vertical diffusion equations¶

7.6.2.4. Solution of the vertical diffusion equations¶

7.6.3. Gravity Wave Drag¶

7.6.3.1. Adiabatic inviscid formulation¶

7.6.3.2. Saturation condition¶

7.6.3.3. Diffusive damping¶

7.6.3.4. Transport due to dissipating waves¶

7.6.3.5. Heating due to wave dissipation¶

7.6.3.6. Orographic source function¶

7.6.3.7. Non-orographic source functions¶

7.6.3.8. Numerical approximations¶

7.6.4. Turbulent Mountain Stress¶

7.6.5. QBO Forcing¶

7.6.6. Radiation¶

7.6.6.1. Longwave radiation¶

7.6.6.2. Shortwave radiation¶

7.6.6.3. Volcanic Heating¶

7.6.7. chemistry¶

7.6.7.1. Chemical Mechanism (Neutral Species)¶

7.6.7.2. Stratospheric Aerosols¶

7.6.7.3. Sedimentation of Stratospheric Aerosols¶

7.6.7.4. Ion Chemistry¶

7.6.8. Electric Field¶

7.6.8.1. Low- and midlatitude electric potential model¶

7.6.8.2. High–latitude electric potential model¶

7.6.8.3. Combing low–/ mid–latitude with the high latitude electric potential¶

7.6.8.4. Calculation of electric field¶

7.6.8.5. Calculation of electrodynamic drift velocity¶

7.6.8.6. Ion drag calculation¶

7.6.9. Boundary Conditions¶