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.2. CAM-Chem using the Modal Aerosol Model

CAM-Chem has the ability to run with two modal aerosols models, the MAM3 and MAM7 ([LEG+12]). The Modal Aerosols Model, is described in Section 4.8.2. In CAM5-Chem, the gase-phase chemistry is coupled to Modal Aerosol Model in chemical species O_{3}, OH, HO_{2} and NO_{3}, as derived from the chemical mechanism and not from a climatoloty. The tropospheric gas-phase and heterogeneous reactions as discribed in Section 4.8.2. are added to the standard MAM chemical mechanism.

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} \cdotF\cdot (1- exp(- \lambda ~ \Deltat)). In this formula, X_{iscav} is the species mass (in kg) and X_{i} scavenged in time step \Deltat , 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:shortwave]).

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_2O), 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
CH_{2}O & x & x &
CO & x & x &
DMS & & x &
ISOP & & & x
NO & x & &
SO_{2} & 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^oN). 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.

  1. update the bottom level values of u, v, q and s using the surface fluxes;
  2. invert (17) and (19) for u,v_{n+1};
  3. compute D and use to update the s profile;
  4. 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 \mum. 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 jth 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 \mum. 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\mum, 10 nm between 5\mum and 10\mum, and 50 nm between 10\mum and 100 \mum.

In the troposphere, stratosphere and lower mesosphere (z < 65km) WACCM4.0 retains the CAM4 shortwave heating (200 nm to 4.55 \mum) 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 jth 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 jth 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 \mum), 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^3K/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\times10^{-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_2O 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_2O reactant (i.e., N_2O_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_2O 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 \mum to approximately 0.5 \mum. 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 \mum; 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_2O is derived in the prognotic water routines of CAM and passed into the chemistry module. Using this condensed phase H_2O, 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\mum. 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_2O 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.

cam6_scientific_guide/figures/ions.jpg

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

Global mean distribution of charged constituents during July solar minimum conditions.
../_images/ionprod1.jpg

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^\circN. 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 \mum, 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(^3P) Implicit   ubvmr
2 O1D O(^1D) 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_2O_2 Implicit dry, wet  
12 N N Implicit   ubvmr
13 N2D N(^2D) 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_2O_5 Implicit    
19 HNO3 HNO_3 Implicit dry, wet  
20 HO2NO2 HO_2NO_2 Implicit dry, wet  
21 CL Cl Implicit    
22 CLO ClO Implicit    
23 CL2 Cl_2 Implicit    
24 OCLO OClO Implicit    
25 CL2O2 Cl_2O_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_3O_2 Implicit    
37 CH3OOH CH_3OOH Implicit dry, wet  
38 CH2O CH_2O Implicit dry, wet flux
39 CO CO Explicit dry flux, ubvmr, airflux
40 CH3CL CH_3Cl Explicit   vmr
41 CH3BR CH_3Br Explicit   vmr
42 CFC11 CFCl_3 Explicit   vmr
43 CFC12 CF_2Cl_2 Explicit   vmr
44 CFC113 CCl_2FCClF_2 Explicit   vmr
45 HCFC22 CHClF_2 Explicit   vmr
46 CCL4 CCl_4 Explicit   vmr
47 CH3CCL3 CH_3CCl_3 Explicit   vmr
48 CF2CLBR CBr_2F_2 (Halon-1211) Explicit   vmr
49 CF3BR CBrF_3 (Halon-1301) Explicit   vmr
50 H2O H_2O Explicit   flux
51 N2O N_2O 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(^1D) + N_2 \rightarrow O + N_2 JPL-06
15 O(^1D)+ O_2 \rightarrow O + O_2(^1\Sigma) JPL-06
16 O(^1D)+ O_2 \rightarrow O + O_2 JPL-06
17 O(^1D)+ H_2O \rightarrow 2 OH JPL-06
18 O(^1D) + N_2O \rightarrow 2 NO JPL-06
19 O(^1D) + N_2O \rightarrow N_2 + O_2 JPL-06
20 O(^1D) + O_3 \rightarrow 2 O_2 JPL-06
21 O(^1D) + CFC11 \rightarrow 3 Cl JPL-06; Bloomfield (1994)
    for quenching of O(^1D)
22 O(^1D) + CFC12 \rightarrow 2 Cl JPL-06; Bloomfield (1994)
23 O(^1D) + CFC113 \rightarrow 3 Cl JPL-06; Bloomfield (1994)
24 O(^1D) + HCFC22 \rightarrow Cl JPL-06; Bloomfield (1994)
25 O(^1D) + CCl_4 \rightarrow 4 Cl JPL-06
26 O(^1D) + CH_3Br \rightarrow Br JPL-06
27 O(^1D) + CF_2ClBr \rightarrow Cl + Br JPL-06
28 O(^1D) + CF_3Br \rightarrow Br JPL-06
29 O(^1D) + CH_4 \rightarrow CH_3O_2 + OH JPL-06
30 O(^1D) + CH_4 \rightarrow CH_2O + H + HO_2 JPL-06
31 O(^1D) + CH_4 \rightarrow CH_2O + H_2 JPL-06
32 O(^1D) + H_2 \rightarrow H + OH JPL-06
33 O(^1D) + HCl \rightarrow Cl + OH JPL-06
34 O(^1D) + HBr \rightarrow Br + OH JPL-06

Table:  Gas-phase Reactions.


no. Reactions Comments
  Nitrogen Radicals  
35 N(^2D) + O_2 \rightarrow NO + O(^1D) \qquad \qquad JPL-06 \qquad \qquad
36 N(^2D) + 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_2O + 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_2O5 + M JPL-06
47 N_2O_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_2O JPL-06
50 NO_2 + HO_2 + M \rightarrow HO_2NO_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_2NO_2 + OH \rightarrow NO_2 + H_2O + O_2 JPL-06
56 HO_2NO_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_2O + 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_2O + O_2 JPL-06
65 OH + OH \rightarrow H_2O + O JPL-06
66 OH + OH + M \rightarrow H_2O_2 + M JPL-06
67 OH + H_2 \rightarrow H_2O + H JPL-06
68 OH + H_2O_2 \rightarrow H_2O + 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_2O_2 + O_2 JPL-06
72 H_2O_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_2O_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_2O \rightarrow HCl + HO_2 + CO JPL-06
79 Cl + CH_4 \rightarrow CH_3O_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_2O_2 + M JPL-06
90 Cl_2O_2 + M \rightarrow 2 ClO + M JPL-06
91 HCl + OH \rightarrow H_2O + 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_2O 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_2O \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_2O 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_3Cl + Cl \rightarrow HO_2 + CO + 2HCl JPL-06
116 CH_3Cl + OH \rightarrow Cl + H_2O + HO_2 JPL-06
117 CH_3CCl_3 + OH \rightarrow 3 Cl + H_2O JPL-06
118 HCFC22 + OH \rightarrow Cl + H_2O + HO_2 JPL-06
119 CH_3Br + OH \rightarrow Br + H_2O + HO_2 JPL-06
  CH_4 and Derivatives  
120 CH_4 + OH \rightarrow CH_3O_2 + H_2O \qquad \qquad JPL-06 \qquad \qquad
121 CH_3O_2 + NO \rightarrow CH_2O + NO_2 + HO_2 JPL-06
122 CH_3O_2 + HO_2 \rightarrow CH_3OOH + O_2 JPL-06
123 CH_3OOH + OH \rightarrow 0.7 CH_3O_2 + 0.3 OH + 0.3 CH_2O + H_2O JPL-06
124 CH_2O + NO_3 \rightarrow CO + HO_2 + HNO_3 JPL-06
125 CH_2O + OH \rightarrow CO + H_2O + H JPL-06
126 CH_2O + 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_2O_5 + H_2O \rightarrow 2 HNO_3 JPL-06; f (sulfuric acid wt %)
2 ClONO_2 + H_2O \rightarrow HOCl + HNO_3 JPL-06; f (T, P, HCl, H_2O, r)
3 BrONO_2 + H_2O \rightarrow HOBr + HNO_3 JPL-06; f (T, P, H_2O, r)
4 ClONO_2 + HCl \rightarrow Cl_2 + HNO_3 JPL-06; f (T, P, HCl, H_2O, r)
5 HOCl + HCl \rightarrow Cl_2 + H_2O JPL-06; f (T, P, HCl, HCl, H_2O, r)
6 HOBr + HCl \rightarrow BrCl + H_2O JPL-06; f (T, P, HCl, HOBr, H_2O, r)
  NAT Aerosol  
7 N_2O_5 + H_2O \rightarrow 2 HNO_3 JPL-06; \gamma = 4 \times 10^{-4}
8 ClONO_2 + H_2O \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_2O JPL-06; \gamma = 0.1
11 BrONO_2 + H_2O \rightarrow HOBr + HNO_3 JPL-06; \gamma = 0.3
  Water-Ice Aerosol  
12 N_2O_5 + H_2O \rightarrow 2 HNO_3 JPL-06; \gamma = 0.02
13 ClONO_2 + H_2O \rightarrow HOCl + HNO_3 JPL-06; \gamma = 0.3
14 BrONO_2 + H_2O \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_2O JPL-06; \gamma = 0.2
17 HOBr + HCl \rightarrow BrCl + H_2O JPL-06; \gamma = 0.3

Table:  Heterogeneous Reactions on liquid and solid aerosols.


no. Reactants Products Comments
1 O_2 + h\nu O + O(^1D) 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(^1D) + 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_2O + h\nu O(^1D) + 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_2O_5 + h\nu NO_2 + NO_3 JPL-06
10 N_2O5 + 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_2NO_2 + h\nu OH + NO_3 JPL-06
15 HO_2NO_2 + h\nu NO_2 + HO_2 JPL-06
16 CH_3OOH + h\nu CH_2O + H + OH JPL-06
17 CH_2O + h\nu CO + 2 H JPL-06
18 CH_2O + h\nu CO + H_2 JPL-06
19 H_2O + 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_2O + h\nu H_2 + O(^1D) (see above)
21 H_2O + h\nu H + 2 O (see above)
22 H_2O_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_2O_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_3Cl + h\nu Cl + CH_3O_2 JPL-06
37 CCl_4 + h\nu 4 Cl JPL-06
38 CH_3CCl3 + 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_3Br + h\nu Br + CH_3O_2 JPL-06
44 CF_3Br + h\nu Br JPL-06
45 CF_2ClBr + 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_3O_2 \sigma: JPL-06;
      based on Brownsword et al. (1997)
48 CH_4 + h\nu H_2 + 0.18 CH_2O + 0.18 O  
   
  • 0.44 CO_2 + 0.44 H_2
see above
   
  • 0.38 CO + 0.05 H_2O
 

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(^2D) 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(^2D) + O 82.389
O_2^+ + e \rightarrow 1.15O + 0.85O(^1D) 508.95
N_2^+ + e \rightarrow 1.1N + 0.9N(^2D) 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(^2D) + N^+ + e
N_2 + e^* \rightarrow N(^2D) + 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_2O, CO_2, H_2, CFC-11, CFC-12, CFC-113, HCFC-22, H-1211, H-1301, CCl_4, CH_3CCH_3, CH_3Cl, and CH_3Br 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_2O, 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.