Coupled Stratospheric Chemistry–Meteorology Data Assimilation. Part I: Physical Background and Coupled Modeling Aspects

Abstract

A coupled stratospheric chemistry–meteorology model was developed by combining the Canadian operational weather prediction model Global Environmental Multiscale (GEM) with a comprehensive stratospheric photochemistry model from the Belgian Assimilation System for Chemical ObsErvations (BASCOE). The coupled model was called GEM-BACH for GEM-Belgian Atmospheric CHemistry. The coupling was made across a chemical interface that preserves time-splitting while being modular, allowing GEM to run with or without chemistry. An evaluation of the coupling was performed by comparing the coupled model, refreshed by meteorological analyses every 6 h, against the standard offline chemical transport model (CTM) approach. Results show that the dynamical meteorological consistency between meteorological analysis times far outweighs the error created by the jump resulting from the meteorological analysis increments at regular time intervals, irrespective of whether a 3D-Var or 4D-Var meteorological analysis is used. Arguments in favor of using the same horizontal resolution for chemistry, meteorology, and meteorological and chemical analysis increments are also presented. GEM-BACH forecasts refreshed by meteorological analyses every 6 h were compared against independent measurements of temperature, long-lived species, ozone and water vapor. The comparison showed a relatively good agreement throughout the stratosphere except for an upper-level warm temperature bias and an ozone deficit of nearly 15%. In particular, the coupled model simulation during an ozone hole event gives better ozone concentrations than a 4D-Var chemical assimilation at a lower resolution.

1. Introduction

This two-part paper describes and examines the benefits of a coupled meteorology–stratospheric chemistry data assimilation system, using the Canadian Meteorological Center (CMC) Numerical Weather Prediction (NWP) model and a comprehensive stratospheric chemistry model, BASCOE (Belgian Assimilation System for Chemical ObsEvations), both validated in their respective environment (i.e., meteorology/chemistry) with data assimilation capabilities. Assimilation experiments, using both weak and strong data assimilation coupling, are presented in Part II [1]. Here in Part I we cover primarily topics relevant to Part II. We examine; 1—the physical basis of stratospheric chemistry–meteorology interactions, 2—model coupling aspects such as the resolution of each geophysical model component (meteorology, chemistry) and whether it should be offline or online in order to yield the most accurate chemical representation, and 3—an evaluation of the coupled model in the context of data assimilation.

The last item deserves some explanation. In model evaluation we should distinguish between unconstrained climate models and atmospheric models where the stratospheric dynamics is constrained by the assimilation of temperature (i.e., radiances). In climate models we compare the climatology of the model against the climatology of the observations. However, data assimilation is based on using the instantaneous observed values at their specific locations. We know that the atmospheric dynamics is chaotic, and its model representation is very sensitive to initial conditions. Atmospheric chemistry models are quite different in that respect; the chemistry is strongly dependent on the meteorology as well as chemical sources and sinks. A free chemistry simulation, without chemical data assimilation but driven by meteorological analyses, can in fact be compared to a reasonable degree of accuracy, to chemical observations at their proper time and location [2]. Therefore, this is the framework we want to use in our chemical modeling evaluation. Furthermore, data assimilation schemes are based on the assumption that differences between observations and model values at their proper time and location (called innovations) are essentially random rather than systematic. In the case where systematic differences are comparable to or exceed the random error standard deviation, the assimilation system is sub-optimal and observations can have little impact (and can even be detrimental, such as introducing spurious spatial structures).

Traditionally, Chemical Transport Models (CTM’s) [3,4] have offered the most complete chemical representation and are often used as benchmarks for chemistry. CTM’s are usually driven offline by meteorological analyses. While the meteorological analyses are usually (and practically) only available at 6 or 12 h time intervals, the meteorology has to be interpolated in time in order to give a dynamical field to drive the chemistry at each model time step. Clearly an offline CTM cannot be used to study the impact of the chemistry on the meteorology.

It is interesting to note the complementary nature of Numerical Weather Prediction (NWP) models with CTMs. Indeed, NWP models solve for; (a) momentum, (b) thermodynamics, (c) conservation of mass and aside from water vapor, (d) use climatological fields as input for chemical composition (in particular ozone and greenhouse gases). In contrast, CTMs; (a) provide a comprehensive representation of chemical composition, (b) solve for the conservation of mass of individual species using chemical reactions and photochemistry, (c) but require as input the momentum (winds), thermodynamics (temperature) and total mass (surface pressure). NWP models routinely use data assimilation, and several CTM’s also have (chemical) data assimilation capabilities. By bringing together these two approaches we can develop a fully coupled chemistry–meteorology model with data assimilation capabilities.

Coupled meteorology-chemistry models provide not only a consistent treatment of the processes shared by meteorology and chemistry but also allow for three-way interactions between physical, chemical and radiation processes [5]. Coupled meteorology-chemistry models are used in several areas such as; (1) in climate simulations with Global Chemistry Circulation Models (GCCMs) [6,7,8] and in climate-chemistry process validation [9,10], (2) for air quality modeling and prediction [11,12], (3) to examine the impact of chemical composition on weather prediction [5,13], and (4) in chemical-meteorology reanalysis such as using weak data assimilation coupling with aerosol-radiation interaction [14,15], or with an ensemble data assimilation but with no interaction and no cross-covariances between species or between species and meteorology [16,17], or using 4D-Var but with chemical increments decoupled from meteorology increments [18]. This wide range of modeling activity is also nicely summarized in a WMO GAW report [19].

The above four classes of coupled models differ somewhat from each other. While GCCM’s consider simulations on multi-decadal to century timescales, air quality models and numerical weather prediction models coupled with air quality focus on time scales of hours to a year. Chemical reanalysis models lie somewhere in between, focusing on periods of a year to a decade. Chemistry climate models and GCCM’s typically have a lower horizontal and vertical resolution than NWP models and generally have limited chemistry composition modeling aimed primarily at simulating GHG (greenhouse gases), aerosols, and aerosol precursors such as those involved in the sulfur cycle (e.g., [20]). When used in process studies (e.g., [10]) or in comparison with observations, these models are usually forced towards meteorological analyses. This may be accomplished either by replacing the model dynamical fields by meteorological analyses (which we call meteorological refresh) or by specified dynamics [7] that consist of a linear relaxation technique to force incrementally the dynamical fields towards linearly interpolated (in time) meteorological analyses (usually called analysis nudging) [21] (not discussed here; see e.g., Froidevaux et al., 2019 [22] for a comparison between the free-running and Specified-Dynamics (SD) configurations of WACCM). In all these cases, except for the GCCM, the radiation feedback on meteorology cannot be examined.

In contrast, air quality models and NWP models coupled with air quality models (classes 2 and 3 above) focus on shorter timescales, have higher resolution both horizontally and vertically and have comprehensive chemistry and aerosol physics and dynamics. Such coupled models allow for the direct assimilation of meteorological observations but are generally used for tropospheric applications [5]. Although temperature and winds have an important effect on the chemical transport and composition in the stratosphere, ozone has an impact on lower stratosphere temperature predictability [23].

In this study, we have constructed such a coupled model starting from the Canadian operational meteorological model [24] GEM (Global Environmental Multiscale model), extended it with relevant physical processes in the stratosphere, and combined it with a comprehensive stratospheric chemical transport model, BASCOE (where as an offline model it has been validated using four-dimensional variational assimilation and ensemble Kalman filtering methods for chemical data assimilation [25,26,27,28]). This coupled model is called GEM-BACH for GEM-Belgian Atmospheric CHemistry model.

The organization of Part I can be summarized as follows. First we present the physics of the coupling between dynamics, radiation and chemical composition in the stratosphere (Section 2). Then, in Section 3, we describe the formulation of the coupled model, what changes were needed, and discuss the modular design of the chemical interface that allows the meteorological model to run with or without chemistry. We then discuss how coupled models and CTM models can be driven by meteorological analyses, discuss their properties and errors, and through a series of experiments we quantify the errors in each formulation. We also discuss the importance of analysis and model resolution on the accuracy of a simulation and contrast it with lower-resolution chemical data assimilation. These considerations suggest that the best approach is to drive a coupled model with meteorological analyses. An evaluation of the modeled stratospheric temperatures and chemical species is then carried out, leading to important considerations for the assimilation component addressed in the second part of this study.

2. Physical Background on Dynamical–Photochemical–Radiation Interactions in the Stratosphere

The material presented in this section is general and describes the diverse interactions from a physical point of view. It applies to both Part I and Part II, and particularly in Part II for both weak and strong data assimilation coupling experiments and the prescription of the error covariances.

From the point of view of energy there are two sources which have a profound impact on the temperature and circulation that characterizes the stratosphere; one is of chemical origin and the other of wave/mechanical origin.

Several tropospheric chemical source species which enter the stratosphere (e.g., O₂, H₂O, N₂O, CFC’s) are photo-dissociated by solar ultra-violet light, producing chemically active (fast reacting) species (e.g., O, OH, NO, ClO, ... ; see

Table 1. List of chemical species in BASCOE.

Source species
Natural	H2O, N2O, CH4, CH3Cl, CH3Br
Anthropogenic	CFC-11 (CFCl3), CFC-12 (CF2Cl2), CFC-113, CFC-114, CFC-115,
	HA-1301 (CBrF3), H-1211 (CBrClF2), HCFC-22 (CHClF2),
	CCl4, CH3CCl3, CHClF2, CHBr3
Short-lived species
Oxygen (Ox)	O3, O(1D), O(3P)
Hydrogen (HOx)	H, OH, HO2, H2O2
Nitrogen (NOx)	N, NO, NO2, NO3
Chlorine (ClOx)	ClOO, OClO, Cl, ClO, ClNO2, HOCl, Cl2O2, Cl2
Bromine (BrOx)	Br,Br2,BrO,BrCl,HOBr
Hydrocarbons (HC)	CH3, CH3O, CH3O2, CH2O, CH3OOH
Long-lived species
	HNO3, HNO4, N2O5, ClONO2, BrONO2
	HBr, HCl, CO, HF, HCO, H2

for a more complete list). That is the case for molecular oxygen, which is the second most abundant atmospheric gas. O₂ is photo-dissociated at an altitude of about 50 km, resulting in atomic oxygen and very rapidly recombines with molecular oxygen giving rise to ozone O₃, and O₃ also recombines with atomic oxygen. These “recombination-type” chemical reactions are exothermic and release so much heat that they transform the vertical stratification of the atmosphere into a deep stable layer from the tropopause up to 50 km, which characterizes the stratosphere. These chemical reactions are known collectively as the Chapman mechanism and are identified with (*) in

Table A1. Gas-phase reactions.

(*) O + O2	→	O3	(*) O + O3	→	2O2	O1D + N2	→	O + N2
O1D + N2	→	N2O	O1D + O2	→	O + O2	O1D + O3	→	2O2
O1D + O3	→	O + O + O2	O1D + H2O	→	2OH	O1D + H2	→	OH + H
O1D + CH4	→	CH2O + H2	O1D + CH4	→	CH3 + OH	O1D + N2O	→	O2 + N2
O1D + N2O	→	NO + NO	O + O	→	O2	ClC4 + O1D	→	4Cl
CFC11 + O1D	→	3Cl + HF	CFC12 + O1D	→	2Cl + 2HF	CFC113 + O1D	→	3Cl + 3HF
CFC114 + O1D	→	2Cl + 4HF	CFC115 + O1D	→	Cl + 5HF	HCFC22 + O1D	→	Cl + 2HF
HA1211 + O1D	→	Br + Cl + 2HF	HA1301 + O1D	→	Br + 3HF	CH3Br + O1D	→	Br
HCFC22 + OH	→	Cl + H2O	CH3Cl + OH	→	HO2 + Cl	CH3Cl + Cl	→	2HCl
CH3CCl3 + OH	→	3Cl + H2O	CH3Br + OH	→	Br + H2O	CHBr3 + OH	→	3Br + H2O
H + O2	→	HO2	H + O3	→	OH + O2	H2 + OH	→	H2O + H
OH + O3	→	HO2 + O2	OH + O	→	O2 + H	OH + OH	→	H2O + O
OH + OH	→	H2O2	HO2 + O	→	OH + O2	HO2 + O3	→	OH + 2O2
H + HO2	→	2OH	H + HO2	→	H2O + O	H + HO2	→	H2 + O2
HO2 + OH	→	H2O + O2	HO2 + HO2	→	H2O2 + O2	H2O2 + OH	→	H2O + HO2
H2O2 + O	→	OH + HO2	H2 + O	→	OH + H	NO + O3	→	NO2 + O2
NO + HO2	→	NO2 + OH	NO2 + O	→	NO + O2	NO2 + O	→	NO3 + O2
NO + O	→	NO3 + O2	NO2 + O3	→	NO3 + O2	NO2 + OH	→	HNO3
NO2 + HO2	→	HNO4	NO3 + O	→	O2 + NO2	NO3 + NO	→	2NO2
NO3 + NO2	→	N2O5	N2O5	→	NO2 + NO3	HNO3 + OH	→	H2O + NO3
HNO4 + OH	→	H2O + NO2 + O2	HNO4	→	HO2 + NO2	NO3 + OH	→	NO2 + HO2
NO3 + HO2	→	NO2 + OH + O2	NO3 + HO2	→	HNO3 + O2	N + NO	→	N2 + O
N + O2	→	N2 + O	NO + O	→	NO2 + O	Cl + O2	→	ClOO
Cl + O3	→	ClO + O2	Cl + H2	→	HCl + H	Cl + CH4	→	HCl + CH3
Cl + CH2O	→	HCl + HCO	Cl + HO2	→	HCl + O2	Cl + HO2	→	OH + ClO
Cl + H2O	→	HCl + HO2	Cl + HOCl	→	Cl2 + OH	Cl + HOCl	→	ClO + HCl
Cl + OClO	→	ClO + ClO	Cl + ClOO	→	Cl2 + O2	Cl + ClOO	→	ClO + ClO
ClO + O	→	Cl + O2	ClO + OH	→	HO2 + Cl	ClO + OH	→	HCl + O2
ClO + HO2	→	O2 + HOCl	ClO + NO	→	NO2 + Cl	ClO + NO2	→	ClONO2
ClO + ClO	→	Cl + OClO	ClO + ClO	→	Cl + ClOO	ClO + ClO	→	Cl2 + O2
ClO + ClO	→	Cl2O2	ClOO	→	Cl + O2	ClO + NO3	→	ClOO + NO2
Cl2O2	→	2ClO	HCl + OH	→	H2O + Cl	HCl + O	→	OH + Cl
OClO + O	→	ClO + O2	OClO + OH	→	HOCl + O2	OClO + NO	→	ClO + NO2
HOCl + O	→	ClO + OH	HOCl + OH	→	H2O + ClO	Cl2 + OH	→	HOCl + Cl
ClONO2 + O	→	ClO + NO3	ClONO2 + OH	→	HOCl + NO3	ClONO2 + Cl	→	Cl2 + NO3
NO2 + Cl	→	ClNO2	NO3 + Cl	→	ClO + NO2	Cl2 + O1D	→	ClO + Cl
HCl + O1D	→	OH + Cl	Cl2O2 + Cl	→	Cl2 + Cl + O2	Br + O3	→	BrO + O2
Br + HO2	→	HBr + O2	Br + CH2O	→	HBr + HCO	Br + OClO	→	BrO + ClO
BrO + O	→	Br + O2	BrO + HO2	→	HOBr + O2	BrO + NO	→	Br + NO2
BrO + NO2	→	BrONO2	BrO + ClO	→	Br + OClO	BrO + ClO	→	Br + ClOO
BrO + ClO	→	BrCl + O2	BrO + BrO	→	2Br + O2	BrO + BrO	→	Br2 + O2
HBr + OH	→	Br + H2O	HBr + O	→	Br + OH	HOBr + O	→	BrO + OH
Br2 + OH	→	HOBr + Br	BrO + OH	→	HO2 + Br	HBr + O1D	→	OH + Br
CO + OH	→	H + CO2	CH4 + OH	→	CH3 + H2O	CH2O + OH	→	HCO + H2O
CH2O + O	→	HCO + OH	HCO + O2	→	CO + HO2	CH3 + O2	→	CH3O2
CH3O + O2	→	CH2O + HO2	CH3O2 + NO	→	CH3O + NO2	CH3O2 + HO2	→	CH3OOH + O2
CH3O2 + NO	→	CH3O + NO2	CH3O2 + HO2	→	CH3OOH + O2	CH3OOH + OH	→	CH3O2 + H2O
CH3OOH + OH	→	CH3O2 + H2O	CH3OOH + OH	→	CH2O + HO2 + OH	CH2O + NO3	→	CO + HO2 + HNO3
CO + O	→	CO2

and

Table A3. Photolysis reactions.

(*) O2 + h$\nu$	→	2O	(*) O3 + h$\nu$	→	O + O2	O3 + h$\nu$	→	O1D + O2
HO2 + h$\nu$	→	OH + O	H2O2 + h$\nu$	→	2OH	NO2 + h$\nu$	→	NO + O
NO3 + h$\nu$	→	NO2 + O	NO3 + h$\nu$	→	NO + O2	N2O5 + h$\nu$	→	NO2 + NO3
HNO3 + h$\nu$	→	OH + NO2	HNO4 + h$\nu$	→	OH + NO3	HNO4 + h$\nu$	→	HO2 + NO2
Cl2 + h$\nu$	→	2Cl	OClO + h$\nu$	→	O + ClO	Cl2O2 + h$\nu$	→	Cl + ClOO
HOCl + h$\nu$	→	OH + Cl	ClONO2 + h$\nu$	→	Cl + NO3	ClONO2 + h$\nu$	→	Cl + NO2 + O
ClNO2 + h$\nu$	→	Cl + NO2	BrCl + h$\nu$	→	Br + Cl	BrO + h$\nu$	→	Br + O
HOBr + h$\nu$	→	Br + OH	BrONO2 + h$\nu$	→	Br + NO3	BrONO2 + h$\nu$	→	BrO + NO2
CH2O + h$\nu$	→	HCO + H	CH2O + h$\nu$	→	CO + H2	CH3OOH + h$\nu$	→	CH3O + OH
ClOO + h$\nu$	→	O + ClO

of Appendix A. The ozone chemistry also involves catalytic loss cycles with the hydrogen (HOx), nitrogen (NOx), and halogen (ClOx, BrOx) families, which generate other constituents and in particular stable (long-lived) molecules called reservoir species (e.g., [29]).

Figure 1 is a diagram of the dynamical-photochemical-radiative interactions in the stratosphere. This major source of heating which is associated with the production of ozone is depicted as a pink arrow in the figure (lower left side). On a global, yearly averaged scale, the heating is nearly counterbalanced by infrared cooling by CO₂, and by O₃ (about half of the effect of CO₂), with a small contribution due to H2O [30] (see Figure S1 in Supplementary Material).

The other important source of energy in the stratosphere is mechanically driven by the drag force due to breaking waves of tropospheric origin. In contrast to the troposphere, where diabatic heating creates vertical motion, it is the mechanically driven circulation in the stratosphere which induces vertical motion of which diabatic heating (in isentropic coordinates) is an outcome and not a cause [31,32]. This source of energy is depicted as the second pink arrow in Figure 1 (upper right side).

The wave drag is explained by the breaking of vertically propagating Rossby and gravity waves in the stratosphere [31]. Planetary-scale Rossby waves, which are forced by orography and land-sea contrasts, can propagate upward only in westerly flow and can reach the stratosphere. Depending on the background mean zonal wind and the wavelength, Rossby waves can become stationary and have growing amplitudes. If, in addition, there is shear flow, tongues of potential vorticity [33,34,35,36] develop, and show up in chemical tracer fields as filamentary structures with cutoff features when the Rossby wave is breaking (e.g., as shown with chemical data assimilation [37]). When the waves break, which occurs at the limit of the diffusion length-scale, they transfer their westward momentum to the zonal flow, thus decelerating it—the so-called wave drag.

Gravity waves have the property that their amplitude increases with height as a result of the decreasing air density. Gravity waves excited by orography reach their critical (breaking) level in the upper troposphere and lower stratosphere [38], whereas those induced by non-stationary waves (such as in frontal systems) break at higher altitudes and play a major role in the middle-atmosphere general circulation [39,40]. The horizontal scale of these gravity waves are much smaller than the typical resolution of global models and thus both their generation and impact must be parameterized.

The wave drag induces a meridional circulation which changes temperatures in isentropic coordinates. Since angular momentum is conserved, the zonal momentum balances the wave drag, and thus implies a negative mechanical forcing that needs to be compensated by a poleward meridional mass flux due to the Coriolis effect. By mass conservation, the meridional mass flux is also linked to the vertical mass flux, creating a meridional circulation with ascent near the tropics and descent near the poles. This is the Brewer–Dobson circulation [31,32,41,42] (see white arrow Figure S2 in Supplementary Material). The vertical motion across isentropes might seem at first perplexing, but a slow persistent vertical motion can move air parcels across isentropes due to the relaxation effect of radiation. Indeed, if an air parcel is displaced downward at a given location, the immediate response is to warm it adiabatically. Then, as the temperature locally rises above the radiative equilibrium temperature, it experiences infrared cooling, which allows the downward displacement to continue. In an isentropic vertical coordinate system (i.e., in a coordinate based on potential temperature) the vertical velocity simply equals the net diabatic heating. The stratospheric meridional circulation that is driven by forces of tropospheric origin pulls the middle atmosphere away from the radiative equilibrium locally, but not globally on long time scales.

2.1. Ozone-Temperature Interaction

Ozone and temperature are related through radiation and photochemistry, but each process results in different ozone-temperature correlations and has its own timescales (see right and left horizontal green arrows in Figure 1). We will discuss first the processes and then their timescales.

The absorption of solar UV radiation in the production of ozone creates a rapid, local increase in temperature, which pulls the temperature away from radiative equilibrium. The perturbed air parcel then undergoes infrared cooling on a slower timescale, ${\tau }_{rad}$, and an adjustment towards a new but higher equilibrium temperature value takes place. Thus, there is a positive correlation between O₃ and temperature because of radiative coupling.

The photochemistry gives a different correlation. Since chemical reaction rates depend on temperature, the ozone production rate increases with decreasing temperature. In terms of absolute value, correlations as high as 0.9 have been reported, based on ozone and temperature measurements from MLS [43] and CRISTA [44]. Several authors have pointed out [43,45,46,47] that the temperature dependence can be represented by a function of the form,

$${\mathrm{O}}_3=Bexp\left(\frac{\Theta }{T}\right)$$

(1)

where B and $\Theta$ are constants whose values depend on the reactants involved in the photochemistry. Taking the derivative of Equation (1) we get the perturbation equation,

$$\frac{\Delta {\mathrm{O}}_3}{{\mathrm{O}}_3}=-\frac{\Theta }{T^2}\Delta T$$

(2)

which shows that temperature perturbations give rise to negatively correlated ozone perturbations.

Let us now discuss the timescales of the different processes. To estimate the radiative timescale ${\tau }_{rad}$ we can use the Newtonian cooling approximation to compute the time required to cool an air parcel through IR emission out to space. This approximation is generally valid above ∼25 km (or 25 hPa) where the radiation exchange between layers can be neglected. For small temperature perturbations, this Cool-to-Space process $Q^{CtS}$ can be written as $Q^{CtS}=-\alpha (T-T_{rad})$ where $T_{rad}$ is a reference temperature near radiative equilibrium and ${\tau }_{rad}=1/\alpha$ is the radiative relaxation timescale. This parameter has been estimated with the GEM-BACH model using a centered finite difference expression ${\tau }_{rad}=2\delta T/[Q(T+\delta T)-Q(T-\delta T)]$ [48] and is displayed in the panel (a) of Figure 2 as a function of latitude and height for a given summer day. Please note that the lower the temperature, the longer is the radiative timescale, thus ${\tau }_{rad}$ decreases with altitude so that a rapid adjustment of the temperature perturbations by the Cool-to-Space process occurs. In the lower stratosphere below (below ∼20 hPa), the radiative timescale is on the order of one month or more, indicating that the temperature perturbations can accumulate over that time period and consequently produce a significant temperature response. The lower stratosphere is thus sensitive to ozone–radiation perturbations.

The ozone photochemical lifetime ${\tau }_{O_3}$ can be defined as the time required for ozone to be reduced by a factor $1/e$ through photochemical reactions. This timescale is readily available from linearized ozone chemistry schemes such as the LINOZ model [49], $d{O}_3/dt=c_1+c_2(O_3-\overline{O_3})+c_3(T-\overline{T})+c_4(O_3^{\uparrow }-\overline{O_3^{\uparrow }})$ which represents the tendency of daily mean values. The overbar denotes the climatology, the ↑ denotes the overhead column of ozone, and the coefficients $c_1,c_2,c_3,c_4$ are determined using a chemical box model. The photochemical timescale for ozone is then ${\tau }_{O_3}=1/c_2$, which is plotted in panel (b) of Figure 2 for 1 June conditions. Please note that the region south of 60° S is in polar night, and the photochemical timescale is infinite.

Comparing the photochemical timescale ${\tau }_{O_3}$ with the radiative timescale ${\tau }_{rad}$ in sunlight conditions we note that above ∼10 hPa, ${\tau }_{O_3}<<{\tau }_{rad}$. The ozone–radiation feedback is small because the lifetime of ozone perturbations is too short to have a significant radiative effect. This is the photochemistry-dominated region, in which ozone and temperature perturbations are negatively correlated. In this region the assumptions for chemical transport modeling are valid, as there is no need to change the temperatures. Below ∼10 hPa, we have ${\tau }_{O_3}>>{\tau }_{rad}$ and the ozone behaves as a passive tracer since the photochemistry can be neglected, but the radiative forcing associated with ozone perturbations persists over several weeks. The ozone radiative impact on temperature can be significant even though the radiative forcing itself is small. Thus, below ∼10 hPa is the radiation-dominated region where ozone and temperature perturbations are positively correlated and there are benefits from a coupled radiation-chemistry model approach.

2.2. Temperature-Wind Interaction

Horizontally, most of the stratosphere is in geostrophic balance, except between about 20° N and 20° S and the upper-stratosphere and mesosphere due to gravity wave breaking. For example, the geostrophic winds derived from satellite observations provided by the CRISTA instrument showed that on a day-to-day basis these winds are remarkably close to the stratospheric winds in the UKMO meteorological analysis [50]. From a dynamical perspective, local perturbations of the horizontal wind and temperature adjust to a balanced state on short time scales (typically less than six hours) by dispersing away fast-moving inertia-gravity waves—a process known as geostrophic adjustment. A scale analysis using the shallow water model reveals that the Rossby number, defined as $R_o=U/fL$, determines the type of adjustment that will take place: when $R_o<1$, the temperature tends to adjust to the wind field, and when $R_o>1$, the wind field tends to adjust to the temperature field (f is the Coriolis parameter and L the length-scale of the disturbance). For planetary-scale waves, such as vertically propagating Rossby waves that enter the stratosphere, $R_o<1$ so that the wind field adjusts to the temperature field. On the other hand, gravity waves of tropospheric origin generally have $R_o>1$, so that the temperature is adjusted to the wind field. Because this adjustment process is mostly completed after six hours, the short term forecast error used in an intermittent assimilation cycle, are in geostrophic balance. In geostropic balanced flow the vertical rate of change of the wind is related to the horizontal temperature gradient by the so-called thermal–wind relation,

$$\frac{\partial u}{\partial p}=\frac{R}{fp}{\left(\frac{\partial T}{\partial y}\right)}_p;\frac{\partial v}{\partial p}=-\frac{R}{fp}{\left(\frac{\partial T}{\partial x}\right)}_p,$$

(3)

where u and v are the zonal and meridional wind components, p is the pressure, R the gas constant, and T the temperature. The thermal–wind relation introduces a three-dimensional coupling between temperature and winds, and the timescale associated with this coupling is on the order of the geostrophic adjustment timescale, ${\tau }_{geos}$ (depicted as the slanted leftward green arrow in Figure 1).

On a much longer timescale and as a result of the Brewer–Dobson circulation, the vertical wind and temperature are also related, but as explained earlier (beginning of Section 2), it is the vertical motion, induced by wave breaking, that determines the temperature distribution. Since the radiation relaxation takes place on a timescale which is faster than the residual circulation, the temperature adapts to a new radiative equilibrium as the fluid particles rise in the tropics or descend in the polar regions. On this slow timescale, we thus observe that the vertical motion drives the temperature change and not the reverse. The timescale ${\tau }_{dyn}$ appearing in Figure 1 refers to the slow timescale associated with the residual circulation.

2.3. Wind-Tracer Interaction

Winds drive the transport of chemical species. From a physical point of view, the chemical tracer mixing ratios have no impact on the winds. For ozone there could be an indirect impact through the radiation followed by geostrophic adjustment, but it is known to be a small, second-order effect, considered to be negligible [51,52].

It has also been observed from aircraft, balloon and satellite platforms that long-lived species display compact relationships in concentrations between species. Also, species with very different sources and sinks exhibit nearly identical meridional-vertical isopleth shapes, indicating that it is the atmospheric transport which maintains these relationships. It was argued [53] that if the sources and sinks are sufficiently slow compared with dynamical timescales, then the meridional slopes of mixing ratio isopleth and the compact correlations between different species are determined by quasi-horizontal mixing. Indeed, the tendency to flatten the isopleths that results from quasi-horizontal mixing is larger than the mean overturning circulation (i.e., the Brewer–Dobson circulation), which tends to steepen the isopleths [53,54]. Depending on the process being considered, the timescales of the tracer-wind relationship occur on a wide range of values; from $U/\Delta x$ if no mixing is considered, to timescale longer than a few weeks, but smaller than the Brewer–Dobson circulation timescale.

3. Description of the Coupled Meteorology-Chemistry Model

An online stratospheric chemistry–meteorology model was developed by combining the Canadian GEM (Global Environmental Multiscale) NWP model with the Belgian stratospheric offline chemical transport model (CTM) [24,25,26,55,56]. The NWP model version used in this study with a lid at 0.1 hPa, has been described and validated by Charron et al. [57] using a variational assimilation system for the assimilation of weather variables [55,56]. This model version became operational for NWP in 2009, but relied on a climatology of stratospheric ozone and did not include any representation of stratospheric chemistry (except for a simple parameterization of water vapor).

The Belgian CTM had a comprehensive stratospheric chemistry and had been delivering operational 4D-Var chemical analyses and forecasts for several years prior to the start of this study. The Belgian operational chemical analysis and forecast system is known as BASCOE and the resulting coupled model was named GEM-BACH (GEM with the Belgian Atmospheric CHemistry). The reader should note, however, that in the years our study took place, only a tropospheric version of GEM (with top at 10 hPa) was operational, so we conducted, ourselves, the tropospheric-stratospheric meteorological assimilation using AMSU-A channels 11–14 for the stratosphere.

Both modeling systems had the necessary components to produce a realistic representation of the stratosphere. The chemistry of the BASCOE CTM, which used a flux-form semi-Lagrangian method for transport [58] was extracted and implemented through a chemical interface in GEM. The resulting chemical transport was then a semi-Lagrangian advection. The stratospheric chemistry component (discussed in Section 3.2 and Section 3.3) and the chemical interface (described in Section 3.4) made the final model fully integrated, from both a dynamical and ozone–radiation interaction perspective, as well as modular (e.g., allowing switching between different chemical packages or running without chemistry).

3.1. Meteorological Model

GEM is a two time-level semi-Lagrangian fully implicit non-hydrostatic grid point model [24] which uses either uniform or variable horizontal-resolution grids [59]. An Arakawa C discretization is used in the horizontal and a hybrid vertical coordinate with non-staggered finite differences is used in the vertical. The model solves for horizontal and vertical momentum, thermodynamics, continuity, and the advection of an arbitrary number of tracers. For the stratosphere, additional physical parameterizations were added including a non-orographic gravity wave drag (GWD) [39,40] scheme for representing the breaking of unresolved non-stationary waves which has a major impact on the middle-atmosphere general circulation. For the computation of radiative transfer the model uses a radiation scheme based on the k-correlated distribution method [60] which is an efficient approach to compute the radiation using a limited number of absorption coefficients. The radiation scheme also uses an ozone climatology based on Fortuin and Kelder (1998) [61] which is merged with HALOE observations above 0.5 hPa. Finally, a stratospheric water vapor parameterized is employed, via a simple methane oxidation mechanism [62,63]. The description and validation of the stratospheric extension of GEM is discussed in Charron et al. (2012) [57].

GEM-BACH used here is configured to run in hydrostatic mode with a global uniform resolution of 1.5° × 1.5°, i.e., 240 × 120 grid points. The vertical domain extends to 0.1 hPa and uses 80 vertical levels, 27 of which are in the stratosphere. Such high vertical resolution is uncommon to most Global Chemistry Climate Models (GCCM’s) which is a distinctive feature of the coupled model GEM-BACH. Figure 3 compares the vertical grid used in this study with the ECMWF vertical grid used back in 2005 (at the beginning of this study). Altitude and vertical grid spacing are estimated using log-pressure altitudes ($z^{*}=Hln(p_0/p)$), where the surface pressure $p_0$ is set to 1000 hPa and the scale height H is set to 7 km.

A series of meteorology-only climate simulations at a resolution of 1.5°× 1.5° using the additional stratospheric parameterizations (outlined above) was performed and the results were compared against either ERA 40 climatology or observations (see Figures S3–S5 in the Supplementary Material). We observed that the native run (a version of the model with a top at 0.1 hPa but with no stratospheric physical parameterizations), exhibit a tropical tropopause and polar winter stratosphere which were too cold while the summer pole stratosphere was too close to radiative equilibrium. But with these stratospheric parameterizations these issues were significantly corrected compared with the ERA 40 reanalysis. Corrections to the zonal winds were also observed with these additional parameterizations (Figure S4 Supplementary Material), where the Hines GWD scheme reduced the mesospheric jets and the new radiation scheme intensified the zonal wind in the stratosphere. The representation of interannual variability in the zonal winds, obtained from a zonal wind time series in the tropics, was also reasonably well captured with these additional parameterizations (Figure S5 in Supplementary Material).

3.2. Stratospheric Chemistry Model

The photochemical module that has been implemented in GEM-BACH is the one used by the Belgian Institute for Space Aeronomy (BIRA-IASB) which includes 57 species, interacting through 143 gas-phase reactions, 48 photolysis reactions and 9 heterogeneous reactions. Table 1 gives the list of species and Appendix A details the list of photochemical reactions. Appendix B explains the physics of the photochemical reaction rates, the so-called J values. Appendix C gives the lower chemical boundary condition at 400 hPa, the level blow which the chemistry solver is not active. The chemical reaction rates and photodissociation rates follow the Jet Propulsion Laboratory compilation [64]. A complete description of the stratospheric chemistry is not intended to be covered here, and we refer the interested reader to appropriate review papers (e.g., [29]).

The rates for gas-phase and heterogeneous chemistry depend on temperature. A reaction between two molecules has a reaction rate $k_g$ of the form

$$k_g=A{e}^{E/RT},$$

(4)

where E represents the energy of activation, R is the gas constant and A the Arrhenius factor. Reactions involving three molecules can also be pressure-dependent and require more complex formulations. Reaction rates involving aerosols are generally expressed as

$$k_{ae}=\frac{\gamma }{4}{\left(\frac{8kT}{\pi M}\right)}^{1/2}{A}_{ae}$$

(5)

where the term in parentheses represents the molecular mean speed of the gas-phase molecules, which depends on temperature T, the molecular mass M, and the Boltzman constant k. $A_{ae}$ is the aerosol surface area per unit volume and $\gamma$ is the reaction efficiency representing the probability that a reaction takes place following the collision of the molecule with the particle. Chemical rate coefficients are determined experimentally and tabulated for different conditions [64].

Heterogeneous chemistry plays an important role, especially in polar regions, and has been explicitly taken into account in GEM-BACH (see

Table A2. Heterogeneous reactions.

ClONO2 + H2O	→	HOCl + HNO3c	ClONO2 + HClc	→	Cl2 + HNO3c	N2O5 + H2O	→	2HNO3c
N2O5 + HClc	→	ClNO2 + HNO3c	HOCl + HCl	→	Cl2 + H2O	BrONO2 + H2O	→	HOBr + HNO3
HOBr + HCl	→	BrCl + H2O	HOBr + HBr	→	Br2 + H2O	BrONO2 + HCl	→	BrCl + HNO3

). Hydrolysis reactions on the surface of Stratospheric Sulfate Aerosols (SSA) contribute mainly to the removal of active nitrogen in the lower stratosphere. In polar regions, another important class of aerosols is Polar Stratospheric Clouds (PSC). Such clouds usually form from SSA particles and grow at cold temperatures from the uptake of water vapor and nitric acid (see [65] for a review). In GEM-BACH, the surface area available for heterogeneous reactions is parameterized in a crude manner. Instead of using a costly detailed microphysical calculation, we used a climatology of SSA surface area densities (see Figure S6 in Supplementary Material). Type II PSC particles (primarily composed of water ice) are set to appear at temperatures below 186°K with a surface area density equal to $5\times {10}^{-9}$ cm${}^2$/cm${}^3$. Between 186°K and 194°K, they are replaced by Type Ia PSC particles (primarily composed of Nitric Acid Trihydrate, NAT) with the same surface area density. The parameterization of PSCs also incorporates the impact of PSC sedimentation on water vapor (dehydration) and gaseous HNO₃ (denitrification). Exponential loss is prescribed for these two species, with characteristic times of 9 days for water vapor (at the gridpoints where type II PSCs are present) and 100 days for nitric acid (at gridpoints where type Ia PSCs are present).

Since the onset of heterogeneous chemistry on PSC depends on temperature, it is important that the meteorological model can reach the threshold temperature required. Figure 4 shows the 15-year average temperature over the South Pole region (defined as the area south of 60° S) as function of height and the day of the year. Of course some years are different from others, and the temperature is not completely uniform in the polar vortex, but Figure 4 does indicate that the GEM model reaches temperatures below 190° C.

3.3. Chemical Solver

The photochemical module and solver used here follow closely those of [66] used at the Belgian Institute for Space Aeronomy (BIRA-IASB). The chemical solver acts on number densities (expressed as molecules-$m^{-3}$), not on mixing ratios as in transport. Ignoring the issue of advection for the moment, the chemical tendency on a model grid is of the form,

$$\frac{\partial c_i}{\partial t}=P_i\left(\mathbf{c}\right)-L_i\left(\mathbf{c}\right)c_i={\Phi }_i\left(\mathbf{c}\right),$$

(6)

where $\mathbf{c}=(c_1,c_2,...,c_J)$ is the vector of number densities for J chemical species, $P\left(\mathbf{c}\right)$ is the production term and $L\left(\mathbf{c}\right)$ is the loss term. Together they define, ${\Phi }_i\left(\mathbf{c}\right)$, the chemical transformation operator for species i. Equation (6) is obtained by adding all the product and loss processes for each species ($i=1,\dots ,J$) from the list of chemical reactions given in Appendix A. This results in J chemical tendency equations, fewer than the total number of chemical reactions N, $(J<N)$. An example of such a procedure is given in Tables 1 and 2 of Yudin and Khattatov (2010) [67]. The photochemical transformation Equation (6) actually forms a set of stiff non-linear equations that span a wide range of chemical timescales. In principle this is best integrated with an implicit time-discretization scheme (e.g., Backward Euler scheme),

$${\mathbf{c}}^{n+1}={\mathbf{c}}^n+\Delta t\mathbf{\Phi }\left({\mathbf{c}}^{n+1}\right)$$

(7)

which ensures computational stability. However, to deal with the non-linearity of $\mathbf{\Phi }$, a linearization around the state at time $t_n$ (up to second order) can be made,

$${\mathbf{c}}^{n+1}={\mathbf{c}}^n+\Delta t(\mathbf{\Phi }\left({\mathbf{c}}^n\right)+\mathbf{J}({\mathbf{c}}^{n+1}-{\mathbf{c}}^n)),$$

(8)

where $\mathbf{J}=\partial \Phi /\partial \mathbf{c}$ is the Jacobian of the chemical production and loss terms, leading to a semi-implicit scheme of the form,

$${\mathbf{c}}^{n+1}={\mathbf{c}}^n+{(\mathbf{I}-\Delta t\mathbf{J})}^{-1}\Delta t\mathbf{\Phi }\left({\mathbf{c}}^n\right),$$

(9)

which, although it is not guaranteed to be stable, is usually stable in practice. The resulting equation is then linear. The actual numerical scheme that solves the chemistry is a Rosenbrock solver of third-order that is a variant and generalization of Equation (9), where the time step is subdivided into several internal time steps h (here 3) (see [68,69] Section 16.6, [70]). This solver is made numerically stable through the specification of the coefficient of the Rosenbrock scheme [71]. The Fortran code needed to apply the Rosenbrock solver for the chemical kinetic equations can be built by the Kinetic PreProcessor (KPP) [72] that also determines the appropriate magnitude of h based on a tolerance factor set as 0.1 in the current version of the model. The chemical solver is applied from the model lid to 400 hPa due to the lack of tropospheric chemistry in the model. In the three bottom layers, species mixing ratios are specified to a set of values taken from the SLIMCAT CTM [4] and are shown in

Table A4. Chemical lower boundary conditions.

N2O = 322 ppbv,    CH4 = 1.76 ppmv,    CH3Cl = 544 pptv,    CH3Br = 10.56 pptv

CFC-11 (CFCl3) = 260 pptv,    CFC-12 (CF2Cl2) = 544 pptv

CFC-113 = 79.333 pptv,    CFC-114 = 4.25 pptv,    CFC-115 = 4.25 pptv

HA-1301 (CBrF3) = 3.3 pptv,    H-1211 (CBrClF2) = 4.62 pptv,    CCl4 = 100 pptv

HCFC-22 (CHClF2) = 170 pptv,    CH3CCl3 = 45.333 pptv,    CHBr3 = 1.1733 pptv

O3 = 20 ppbv,    O(1D) = 1.E-21,    O(3P) = 2.E-17

H = 2.E-22,    OH = 1.E-15,    HO2 = 1 pptv,    H2O2 = 2 ppbv,    H2 = 1.E-21

N = 1.E-21,    NO = 1.E-13,    NO2 = 2 pptv,    NO3 = 3.E-14

ClOO = 1.E-21,    OClO = 4.E-15,    Cl = 9.E-19,    ClO = 4.E-14

ClNO2 = 1.E-21,    HOCl = 2.E-13,    Cl2O2 = 6.E-21,    Cl2 = 1.E-21

Br = 3.E-18,    Br2 = 1.E-21,    BrO = 7.E-16,    BrCl = 3.E-16,    HOBr = 3.E-15

CH3 = 1.E-21,    CH3O = 1.E-21,    CH3O2 = 1.E-21,    CH2O = 1.E-21

CH3OOH = 0.649 ppbv

HNO3 = 2 pptv,    HNO4 = 3.E-14,    N2O5 = 2.E-14,    ClONO2 = 1.E-13

BrONO2 = 5.E-16

HBr = 4.E-15,    HCl = 1 pptv,    CO = 15 ppmv,    HF = 1.E-21,    HCO = 1.E-21

CO2 = 380 ppmv

in Appendix C. Species vertical fluxes are null at the model lid.

The execution of the chemistry solver with a semi-Lagrangian transport scheme proceeds as follows. First, a semi-Lagrangian advection is performed on all species by interpolating from the upstream (i.e., departure) points to compute the mixing ratio at the arrival point on the model grid. Please note that all species have the same upstream point and interpolation weights, which are calculated only once. This represents a significant computational savings for the chemical transport. The species mixing ratio is converted into number density and the photochemical tendency for each species is computed on each model grid point. Once the number densities are updated, they are transformed back into mixing ratios for another transport time step or for a call to the physics scheme. The transformation from number density to mixing ratio follows the expression $\chi =(R_{*}T/N_{A}p)c$, where T is the temperature, p the pressure, $R_{*}$ the universal gas constant, and $N_A$ the Avogadro number.

3.4. Model Coupling and Interface

Models are composed of several processes which are integrated either sequentially or in parallel (simultaneously). In sequential processing, for a given time step, the model state is updated after each process and provides input to the next process, until all processes are integrated. Sequential processing is also called time-splitting. In parallel processing, the tendencies of each process are computed simultaneously using the same initial model state. The updated state is computed from the sum of tendencies. Parallel processing is also called process splitting.

Parallel processing is appealing because of its simplicity and ease of implementation for coupled models. However, it has the disadvantage that the stationary solution of the time-discrete equation does not match the stationary solution of the time-continuous equation [73,74]. Sequential processing does not have this problem—it has the same stationary solution as the time-continuous equation. However, the transient solution depends on the order in which the different processes are integrated. The total error is minimized when the processes are ordered from the slowest process to the fastest [75,76].

The meteorological model GEM uses sequential processing and we have followed closely the same approach for the coupled meteorology-chemistry model configuration. In addition, we have adopted a modular design such that the chemistry component can be present (or not) through a chemical interface. This flexibility allows having a meteorology-only or meteorology-chemistry model configuration. However, this flexibility entails a small additional computational cost and maintenance, since some physics routines need to be duplicated and present in the chemical module. To make this clear let us begin by discussing what sequential processing would look like if chemistry were completely integrated with the physics module.

The processes in the meteorological model GEM are updated in the following sequence: 1—Radiation, 2—Advection, 3—Dynamics terms of meteorological variables using a semi-implicit scheme, 4—Surface fluxes and gravity wave drag (orographic and non-orographic), 5—Boundary layer processes and vertical diffusion, 6—Shallow convection, 7—Deep convection, and 8—Microphysics. A coupled meteorology-chemistry model has a 9th process—Chemistry, which involves very fast process that in principle should be solved implicitly but in practice we have chosen to use a semi-implicit approach using a Rosenbrock solver (see Section 3.3). For the purpose of this discussion let us consider only those processes that involve both meteorological and chemical variables. For stratospheric chemistry those are: 1—Radiation, 2- Advection, 5—Vertical diffusion, and 9—Chemistry.

Let $\mathbf{X}$ represent the coupled (augmented) state vector, i.e.,

$$\mathbf{X}=\left(\begin{array}{c}\mathit{\mu }\\ \mathit{\chi }\end{array}\right)$$

(10)

where $\mathit{\mu }$ is the meteorological state vector and $\mathit{\chi }$ the chemical state vector. Then the evolution that matters for the coupled state vector takes the form

$$\begin{array}{c}\hfill \frac{D\mathbf{X}}{Dt}=\mathbf{R}\left(\mathbf{X}\right)+\mathbf{D}\left(\mathbf{X}\right)+\mathbf{\Phi }\left(\mathbf{X}\right)\end{array}$$

(11)

where $D\mathbf{X}/Dt$ represents the material derivative, $\mathbf{R}$ the radiation, $\mathbf{D}$ the vertical diffusion and $\mathbf{\Phi }$ chemical processes.

The radiation $\mathbf{R}\left(\mathbf{X}\right)$ can be either offline or online with the prognostic chemical variables–in particular O₃. In the offline mode, greenhouse gases and a zonal-mean climatology of O₃ are given as input to the radiation. Our ozone climatology is based on Fortuin and Kelder [61] and HALOE observations above 1 mb. Thus, the radiation process takes the form

$$\begin{array}{c}\hfill \mathbf{R}\left(\mathbf{X}\right)=\left(\begin{array}{c}\mathbf{R}({\mathit{\mu }}^n,\overline{{\mathit{\chi }}^c})\\ \mathbf{0}\end{array}\right).\end{array}$$

(12)

Since radiation in the troposphere depends on cloud parameters that are diagnostic, and thus not advected, and also depends on temperature that is advected, it is desirable to compute radiation before advection in order to avoid any mismatch in the fields required as input. The radiation update thus operates on the initial state ${\mathbf{X}}^n$ to create an intermediate state ${\mathbf{X}}_R^{*}$ of the form

$${\mathbf{X}}_R^{*}=\left(\begin{array}{c}{\mathit{\mu }}_R^{*}\\ {\mathit{\chi }}_R^{*}\end{array}\right)={\mathbf{X}}^n+\Delta t\mathbf{R}\left({\mathbf{X}}^n\right)=\left(\begin{array}{c}{\mathit{\mu }}^n+\Delta t\mathbf{R}({\mathit{\mu }}^n,\overline{{\mathit{\chi }}^c})\hfill \\ {\mathit{\chi }}^n\hfill \end{array}\right)$$

(13)

In a fully coupled ozone–radiation configuration that we will consider in Part II of this study [1], the radiation process then takes the form $\mathbf{R}=\mathbf{R}({\mathit{\mu }}^n,{\mathit{\chi }}_{\mathrm{O}3}^n)$.

After radiation, advection is processed on both meteorological and chemical variables using ${\mathbf{X}}_R^{*}$ as the initial state. Without loss of generality, in a semi-Lagrangian scheme we can write the advection update as,

$${\mathbf{X}}_A^{*}=\left(\begin{array}{c}{\mathit{\mu }}_A^{*}\\ {\mathit{\chi }}_A^{*}\end{array}\right)=\left(\begin{array}{c}{\mathit{\mu }}_R^{*}(\mathbf{x}-\mathit{\alpha }\Delta t)\\ {\mathit{\chi }}_R^{*}(\mathbf{x}-\mathit{\alpha }\Delta t)\end{array}\right),$$

(14)

where $\mathbf{x}$ is the spatial coordinate and $\mathit{\alpha }$ the upstream displacement along the trajectory. The next process to consider is the vertical diffusion $\mathbf{D}\left(\mathbf{X}\right)$ that is applied on both meteorological and chemical state variables. The resulting update has the form

$$\begin{array}{c}\hfill {\mathbf{X}}_D^{*}=\left(\begin{array}{c}{\mathit{\mu }}_D^{*}\\ {\mathit{\chi }}_D^{*}\end{array}\right)=\left(\begin{array}{c}{\mathit{\mu }}_A^{*}+\Delta t{\mathbf{D}}_K\left({\mathit{\mu }}_A^{*}\right)\\ {\mathit{\chi }}_A^{*}+\Delta t{\mathbf{D}}_K\left({\mathit{\chi }}_A^{*}\right)\end{array}\right),\end{array}$$

(15)

using diffusion coefficients K computed from meteorological fields that are common to both meteorological and chemical variables. Finally, the coupled state is updated for the chemical processes,

$$\begin{array}{c}\hfill {\mathbf{X}}^{n+1}={\mathbf{X}}_{\Phi }^{*}=\left(\begin{array}{c}{\mathit{\mu }}_{\Phi }^{*}\\ {\mathit{\chi }}_{\Phi }^{*}\end{array}\right)=\left(\begin{array}{c}{\mathit{\mu }}_D^{*}\hfill \\ {\mathit{\chi }}_D^{*}+{(\mathbf{I}-\Delta t{\tilde{\mathbf{J}}}_{\mu })}^{-1}\Delta t{\tilde{\mathbf{\Phi }}}_{\mu }\left({\mathit{\chi }}_D^{*}\right)\hfill \end{array}\right),\end{array}$$

(16)

with a resulting state ${\mathbf{X}}^{n+1}$. Here in Equation (16), the dependence on $\mu$ is actually $(T,p)$, temperature and pressure. We also use the tilde $\widetilde{\left(\right)}$ to emphasize that the Jacobian and chemistry production and loss terms are evaluated in terms of mixing ratio and not in terms of the number density as in Equation (9).

Figure 5 displays the sequence of updates of the coupled (augmented) model state and which information is passed from meteorological to chemical modules (up and down arrows). The equal sign indicates that there are no changes and the horizontal arrows indicate changes due to a specific process update. To simplify, let us first discuss the case where there is no ozone–radiation coupling, i.e., let us ignore for now the dotted upward arrow.

First, infrared radiation does not change the chemical concentrations but does change the temperature. Then, the advection of chemical (prognostic) variables requires information about the displacement of the upstream point, $\mathit{\alpha }$, and the interpolation weights, $\mathbf{w}$, that are computed from the wind and the position of the upstream point. Next, applying the vertical diffusion on chemical variables requires sharing the diffusion coefficients, K, that are computed from the meteorology. Finally, chemistry requires information about temperature and pressure, but does not change the meteorological variables. When this last update is completed, the state at time $t^{n+1}$ is produced.

In a modular implementation, all processes that involve changes in chemical concentrations (in practice, with the exception of advection) can be included in a chemical module. The computation of the changes in chemical concentrations requires exchange of information through a chemical interface. For example, if we duplicate and include the vertical diffusion routine in the chemical module and pass the K coefficients to the chemistry module, the computation of the vertical diffusion of the chemical variables can be done in the chemistry module. Likewise, in principle, the advection of the chemical variables could be performed in the chemistry module if we duplicate the appropriate routines, but in practice it is easier (in terms of code maintenance), to simply pass the chemical concentrations to the advection routine that computes advection to all meteorological and chemical variables at once.

So far, the exchange of information is performed only one way, from meteorology to chemistry. It is then possible and easy to have a meteorology-only configuration separate from a meteorology-chemistry configuration, by simply allowing advection to be performed on an arbitrary number of variables.

The modularity of this approach can be preserved with sequential processing in the case of ozone–radiation interaction. Indeed, after the whole sequence of processes from advection to increasingly faster processes is completed, the prognostic ozone can be passed as the initial condition to the radiation scheme for the next model time step integration (see dotted upward arrow in Figure 5).

Finally, in terms of computational resources, GEM-BACH is about five times slower than GEM (with no chemistry) on a uniform resolution 1.5° × 1.5°, i.e., 240 × 120 grid points, with 80 vertical levels, and running on 16 CPU (MPI 4 nodes, OpenMP 4 CPU). GEM-BACH transport (advection of 57 species) accounts for about 1/4 of the CPU time, the computation of the J-values for about 1/4 of the CPU time, and the Rosenbrock chemical solver about 1/2 of the CPU time.

4. Coupling with Meteorological Analyses

In data assimilation we pay a lot of attention to errors in the model and the observations. In a coupled model, it may be important to examine the coupling strategy and the resolution of the different model components, in order to minimize the size and/or number of errors. In the next three sections we will examine different source of errors in dynamical model coupling with meteorological analyses and will assess their impact. This will guide us in choosing the way forward for the comparison of the chemistry simulations with observations (Section 7) and to conduct the combined chemical and meteorological assimilation in Part II [1].

By way of introduction of the relevant topics, let us examine the time continuity. The use of (time) sequential meteorological analyses, such as in 3D-Var, creates a temporal discontinuity in the meteorological fields. The question arises then as to “What is the implication for the chemical transported fields”.

Depending on the type of chemical model coupling, the implications are different. There exist essentially two methods for dynamical coupling. One is the (offline) chemical transport model (CTM) approach which bypasses the meteorological time discontinuity by using an temporal interpolation between meteorological analyses, and the other is initialization of the meteorological forecast with the analysis, and then allowing this forecast to drive the chemical model. This is done with dynamically coupled meteorology-chemistry models. Specifically:

(CTM)

With an offline chemical transport model (CTM). Meteorological analyses which are usually produced at regular time intervals, e.g., 6 h, can be linearly interpolated in time to drive an offline CTM at each time step. This is usually done by interpolating the horizontal wind, temperature and surface pressure and diagnosing the vertical motion from the divergence of the horizontal wind. Alternatively, the vertical motion can also be computed from the diabatic heating rate, giving the vertical motion in isentropic coordinates [77,78], or

(MR)

With a coupled dynamical meteorology-chemistry model using meteorological refresh (MR). In this case, the coupling is achieved by a direct insertion of meteorological analyses at analysis time. As new analyses are available (e.g., each 6 h), the meteorological variables of the coupled model are reset to the given meteorological analysis values, but in between, the internal dynamics of the coupled model come from the meteorological driver.

Please note that there is also a variant to these approaches, called specified dynamics [7,21] (discussed in the Introduction) that combines the linear interpolation between analyses to force incrementally the dynamical fields of a dynamically coupled meteorology-chemistry model, and thus inherit the properties of both approaches.

CTM coupling has the advantage that there is no discontinuity of the meteorological variables at the analysis times. There is a smooth transition of the meteorology from one meteorological analysis time to the next one. However, because of interpolation in time it has the disadvantage that meteorology is not dynamically consistent between analyses.

The MR mode has the advantage and disadvantage interchanged compared to those of the CTM mode. It has the advantage that during the model time integration between analyses, the meteorology is dynamically consistent. But it has the disadvantage that at the analysis time, there is a discontinuity in the meteorological fields, where the jump is a result of the (meteorological) analysis increment.

Table 2. Effective use of meteorology in CTM and MR modes.

	CTM	MR
At analysis times	Continuous meteorology	Discontinuous meteorology
Between analysis times	Dynamically inconsistent	Dynamically consistent

summarizes how meteorology is effectively used in CTM and MR modes.

For either CTM or MR modes, it is important to note that since there is no chemical data assimilation, the chemical concentrations are time-continuous (both at the meteorological analysis time and between analyses). The absence of changes of concentrations immediately before and immediately after a meteorological analysis time $t_A$, can be written as

$$\chi \left(t_A^{-}\right)=\chi \left(t_A^{+}\right).$$

(17)

However, the time derivative of the concentration at the analysis time $t_A$ is given by

$$\frac{d\chi }{dt}{|}_{t_A^{+}}-\frac{d\chi }{dt}{|}_{t_A^{-}}=[\mathbf{V}\left(t_A^{+}\right)-\mathbf{V}\left(t_A^{-}\right)]·\nabla \chi \left(t_A\right)=\left\{\begin{array}{c}0\mathrm{CTM}\mathrm{mode}\hfill \\ \Delta {\mathbf{V}}^A·\nabla \chi \left(t_A\right)\mathrm{MR}\mathrm{mode}\hfill \end{array}\right.$$

(18)

where $\Delta {\mathbf{V}}^A$ is the wind (or meteorological) analysis increment (see Appendix D for a derivation). The time derivative of the concentration is continuous in CTM mode but discontinuous in MR mode.

Let us now outline some properties of the concentration error. In principle, the accumulation of concentration error can be decomposed into two parts: 1—the accumulation of error between the (meteorological) analyses, and 2—an error at the (meteorological) analysis time. However, since the concentration is time-continuous at analysis time (Equation (A9)), and since the true concentrations should also be continuous, we conclude that the concentration error is continuous at the meteorological analysis time (and of course also between meteorological analyses). Thus, in both CTM and MR modes there is no jump in the concentration error at the (meteorological) analysis time. Also, since the meteorological analysis increment is the same in both CTM and MR experiments, the concentration correction (which can be viewed as an unobserved variable in Equation (5), Part II [1]) should also be the same, provided that the error statistics are the same. We thus conclude that any change in concentration error between CTM or MR modes depends essentially on the error in transport between the meteorological analysis times.

We will present in the following section a comparison of experiments using the CTM and MR modes. However, we can already speculate that the use of linearly interpolated winds (as in CTM mode) would create larger transport (concentration) errors than would a dynamically evolving wind field from a coupled meteorology-chemistry model. Thus, we anticipate that concentration errors in MR mode will be smaller than in a CTM mode. The impact of discontinuity in meteorology and dynamical inconsistency is investigated numerically in Section 5.

5. Comparison between CTM and Coupled Meteorology-Chemistry Model

The distribution of atmospheric constituents is strongly driven by photochemistry, emissions and meteorology. For this reason, chemical models driven by meteorological analyses can, to a certain degree, simulate observed concentrations at their proper time and location. In terms of chemistry such simulations are free model runs, while the constraints arising from observations come only from the meteorological analysis.

Meteorological analyses extending to the stratosphere were obtained for the first time under this study (at the time the experiments were carried out, the operational model GEM had a lid at 10hPa only). They were obtained from the CMC 3D-Var-FGAT [55] (i.e., 3D-Var with First Guess at Appropriate Time) and 4D-Var [56] assimilation systems using the tropospheric meteorological observations used operationally at CMC (see [55] Appendix B) and AMSU-A channels 11-14 for the stratosphere using the same bias correction scheme as for tropospheric radiance observations [79] (a procedure also described in Section 4 of [80]). Please note that an improved bias correction will be presented in Part II [1] Section 4.

Several experiments were carried out to quantify the different sources of errors between GEM-BACH in MR mode and the BASCOE CTM. All MR experiments were performed at 1.5° × 1.5° resolution, i.e., 240 × 120 grid points, with 80 levels (27 levels are in the stratosphere) and a 45 min timestep as described in Section 3.1. In the case of the CTM experiments the same horizontal and vertical resolution were used but with a timestep of 15 min because its flux-form semi-Lagrangian requires satisfying the Courant-Friedrichs-Lewy condition in the meridional direction. All experiments were evaluated over a period of 12 days in late summer 2003.

Figure 6 shows the results of an ozone simulation for CTM and MR modes using different meteorological analyses. The comparison is made against limb sounding ozone observations from the Envisat/MIPAS instrument. The solid lines depict the mean difference, and the dashed curves the standard deviation. The dashed curves are plotted symmetrically with respect to the zero error (vertical solid black line), simply to illustrate the range of $\pm \sigma$ random errors. Also note that all the errors are normalized by the observed values, so that the errors are expressed as a percentage. The differences are computed from interpolating the model values to the observation location, at the observation time of validity. We emphasize that the chemistry model is the same in all experiments, and the results differ only due to the meteorological analysis and its coupling to the chemistry transport. The red and grey curves are results from using MR (meteorological refresh) mode using the CMC 3D-Var-FGAT and 4D-Var meteorological analyses respectively, while the green curves are results based on the CTM mode using 3D-Var-FGAT analyses.

First, we draw the reader’s attention to the error standard deviation (i.e., dashed curves). We note that between 100 hPa and 20 hPa where the photochemical timescale is of the order of several months (i.e., transport is dominant), the standard deviation is significantly larger with the CTM mode and reaches twice the MR values at 100 hPa. The difference between red and green curves is only due to the mode of coupling, CTM vs MR. Also, there is little difference between using 3D-Var-FGAT and 4D-Var meteorological analyses. The main difference arises from whether the MR mode or CTM mode is used. This indicates that the dynamical inconsistency in the time integration window of 6 h is the main source of random error, while using different analyses is of secondary importance.

Similar conclusions can be drawn for the systematic error (solid curves), where we observe that the main difference in error arises from using either MR or CTM mode of coupling rather than using different meteorological analyses (either 3D- or 4D-Var).

Above 20 hPa, where the photochemical timescale is shorter and transport plays a negligible role with respect to photochemistry, there is no difference in error standard deviation but only a slight difference in systematic error, probably due to differences in mean temperatures between the different meteorological analyses.

These results clearly indicate that despite the jump in the wind field at the analysis refresh time, the integration consistency which arises in coupling with MR gives a superior chemistry simulation compared to the time interpolation error and associated dynamical inconsistency in CTM coupling. We conclude that for chemical data assimilation, model error due transport is smaller for coupled models in MR mode compared to offline models in CTM mode.

6. Discussion on the Use of Lower Resolution

It has been argued that the tracer distribution is primarily controlled by the large-scale low-frequency component of the flow [31] and thus, a lower-resolution meteorology could be used in a chemical transport model, and yet produce realistic simulations. A similar argument exists for meteorological data assimilation, where lower-resolution analysis increments produced in a 4D-Var scheme (i.e., incremental 4D-Var) are able to create realistic small-scale structures over time through the atmospheric model [81,82]. The arguments presented for both tracer-wind transport and the incremental 4D-Var are based on models that are dominated by an enstrophy cascade with a $k^{-3}$ spectrum, where k is the wavenumber. Large-scale lower-resolution atmospheric and barotropic vorticity models exhibit such behavior—the sources of energy injected at small wavenumbers cascade down to higher wavenumbers. Nevertheless, in the stratosphere, as simulated by high resolution meteorological models, energy spectra evolve from a steep spectrum, $\propto k^{-3}$, to a shallow spectrum as the height increases [83] resulting from an inverse energy cascade $\propto k^{-5/3}$ at larger wave numbers. We thus expect that GEM-BACH has a similar behavior.

The chemical tracer field has what is called a scalar variance spectrum that has in theory a slope lying between −1 for enstrophy-cascade dynamics to −5/3 for inverse energy cascade dynamics [84]. However, because of mixing barriers and trapping by persistent vortices, it has been argued that the scalar spectral slope can be as steep as −2 in those cases [85]. However, in general, stratospheric observations indicate that slopes of −5/3 in the scalar variance are usually obtained [86,87].

In the context of tracer-wind transport, the argument that low-resolution winds can reproduce the small-scale structure of tracer fields has been challenged by Bartello [84]. He pointed out that lower-resolution models dominated by enstrophy-cascade dynamics can reproduce accurately some fine structure of the tracer field [84,88] using a relatively coarse wind field. However, higher-resolution models, with an inverse energy cascade at smaller scales, create fictitious small-scale structures in the tracer field when coarse resolution winds are used [84].

These results are of direct relevance to coupled stratospheric chemistry–meteorology modeling and data assimilation. In the stratosphere, where the inverse energy cascade is important, the use of lower-resolution analysis increments, e.g., incremental 4D-Var, is expected to result in a loss of information in the stratospheric meteorological analysis. Furthermore, chemical tracer fields driven by lower-resolution analysis increments would also lose accuracy in small-scale structures. For these reasons, we argue that for stratospheric models such as GEM-BACH, we should generate analysis increments at the same resolution as the meteorological model, and also drive the chemistry at the same resolution as the meteorological model, so as not to introduce small-scale errors in both meteorology and chemical tracer fields.

Bartello [84] also noted that the temporal resolution of the advecting velocity would have a timescale $\propto k^{-2/3}$ and thus an increase in spatial resolution requires a corresponding increase in temporal sampling. Linear time interpolation used in CTM’s would then be detrimental, which has been clearly identified in Section 5.

To illustrate our discussion regarding resolution, we have chosen an ozone depletion event and conducted a few experiments over the period August–October 2003. First, we have conducted a 3D-Var-FGAT meteorological assimilation with GEM at 1.5° × 1.5° resolution producing analysis increments at the same resolution, and driving the coupled model GEM-BACH in the MR mode also at the same resolution. The choice of the horizontal resolution (i.e., 1.5° × 1.5°) is guided by the ability to accurately model the ozone hole. We remark that many successful studies of polar ozone depletion use CTM or specified dynamics GCCM with latitude–longitude resolutions of approximately 2° × 2.5° [22]. This resolution became typical for GCCM aiming at modeling the ozone hole as the minimum resolution allowing a realistic representation of the tropical and high-latitude mixing barriers (i.e., polar vortex edges) [89]. It should also be noted that several modeling groups use even coarser resolutions to model polar ozone depletion, e.g., with a 3-D CTM at 2.5° × 3.75° [90] or with a chemistry climate model in nudged (SD) mode at T42 (i.e., approximately 2.8° × 2.8°) [91]. Hence the choice of 1.5° × 1.5° was deemed sufficient to accurately model the ozone hole with GEM-BACH.

Figure 7 panel (a) shows the total column ozone measurements from the TOMS instrument on 30 September 2003. Typical to such events, we note there is a wide range of total ozone amounts (and thus of ozone concentration values) and sharp gradients along the vortex edge. The results of pure chemical simulation with MR and at the resolution 1.5° × 1.5° are presented in panel (b) of Figure 7. We observe a remarkably accurate simulation of the ozone depletion event, with accurate vortex values, sharp gradients along the vortex edge and reasonably well-reproduced mid-latitude surf-zone values. Tropical values (see Figure S12 in Supplementary Material) are lower than observed (we know from Section 7.2.3 that the GEM-BACH model has an ozone deficit problem). Nevertheless, this represents a major accomplishment, taking into account the fact that there is no chemical assimilation in this run.

We now compare the BASCOE CTM (panel c) and 4D-Var chemical data assimilation (panel d), as is done operationally at BIRA, using the operational ECMWF meteorological analysis. Since 4D-Var chemical assimilation is costly, the model and assimilation are performed at a lower resolution, in the case here both at 3.75° × 5°. As in GEM-BACH, heterogeneous chemistry is simulated using prescribed climatological SAD values (see Section 3.2). The CTM simulation with the coarser-grained model overestimates both the vortex and mid-latitude surf-zone values and has weak gradients at the vortex edge. Thus, it appears that lower-resolution winds driving a chemical model at the same resolution cannot reproduce the sharp gradients (i.e., small-scale structures in the chemical field) at the vortex edge. Although this result is simply an illustration, it does not support the claim that the large-scale, low-frequency component of the flow controls the tracer distribution. It is important to note that both the BASCOE CTM and GEM-BACH have no horizontal diffusion that could smooth horizontal gradients. Any horizontal gradient or lack thereof is a result of the driving winds and underlying cascade regime.

Next, 4D-Var assimilation of ozone observations from MIPAS was conducted with a low-resolution BASCOE 4D-Var assimilation system and the result is presented in panel (d). We note lower ozone values in the polar vortex than with the BASCOE CTM at the same low resolution (panel c), but the horizontal extent of these low values (column < 125 DU) is not as large as in TOMS observations (panel a). The extent of the depleted area is actually better simulated by the higher-resolution GEM-BACH model (panel b) then by this MIPAS-based analysis (panel d), where we also observe a weakening of the concentrations outside the polar vortex. This is apparently an effect of the error covariances and the impact of observations near the vortex edge that tend to mix values in and out of the vortex, which can be partly alleviated if a more appropriate covariance model is used [92].

These results argue in favor of using the same resolution for both meteorological and chemical components, as well as analysis increments. Combined with the selection of the Meteorological Refresh mode (see previous Section 5), this provides the justification for the data assimilation setup used in the next Section 7 and in Part II [1].

7. Evaluation of GEM-BACH Driven in Meteorological Refresh (MR) Mode

7.1. MIPAS and HALOE Measurements

The MIPAS instrument on-board the Envisat satellite [93] is a limb sounder which uses a Fourier-transform spectrometer for the detection of emission spectra in the middle and upper atmosphere [93,94,95]. It observes a wide spectral interval throughout the mid-infrared with high spectral resolution, which permits retrievals of pressure in addition to temperature and volume mixing ratio (VMR) of different gases. The instrument provides about 1000 profiles per day (day and night) with a global spatial coverage. The operational ESA retrievals v4.61 that we used here do not use any a priori and thus can be considered to be “pure observations”. The typical root mean square error (RMSE) observation error of MIPAS is about 2° K for temperature [96], 10% VMR error for O₃ [97], 20% VMR error for H₂O, CH₄ and N₂O [98,99]. Other species such HNO₃ and NO₂ have a large relative error that varies considerably with altitude (although with a minimum error of 10% at 22 km for HNO₃ and at 40 km for NO₂) [100,101] and were used in this comparison.

The HALOE instrument on-board the UARS satellite [102,103,104] is a solar occultation instrument which employs a broadband radiometer, and a gas correlation technique specifically to infer aerosol extinction [105]. Each HALOE radiometric profile is divided by the exo-atmospheric signal thus giving a direct measurement of the atmospheric transmission for each channel. HALOE retrieves temperature and mixing ratio of several species. The horizontal coverage is limited to two latitudes on a given day, one at sunrise and the other at sunset. These sun occultation latitudes change gradually over a period of 45 days, as the UARS satellite is on an inclined orbit and undergoes precession. The whole latitudinal coverage is quite complicated but ranges from ∼45° in one hemisphere to ∼80° in the other hemisphere. Every 45 days the hemispheric coverage is inverted through a yaw maneuver of the satellite. The HALOE retrieval of atmospheric constituents is a modified “onion peel” algorithm with no a priori information. The main source of error arises from the absence of pressure measurements. Because of the pointing uncertainty, there is a need to perform a registration of profiles with altitude and pressure, which is done by referencing a meteorological analysis that contains errors [104]. HALOE retrieval V19 used in this study is considered an excellent verification dataset with errors less than 10% for ozone, temperature and H₂O [105,106,107] and 15% for CH₄ [108].

7.2. GEM-BACH Evaluation against Satellite Observations

Data assimilation schemes are based on the assumption that errors (model and observations) are random rather than systematic. These (first) evaluations were performed, in part, to establish the readiness in assimilating MIPAS observations, by assessing whether the differences between the model and observations is dominated by random errors rather than biases.

7.2.1. Temperature

Temperature is important for the chemistry and for the thermal–wind component of the transport. Figure 8 shows the difference between temperature observations and GEM-BACH driven in MR mode (positive differences correspond to observed values being larger than modeled values). Results using MIPAS are displayed in the left panel and HALOE in the right panel. The statistics (mean and standard deviation) use model values interpolated to the proper time and location of the individual observations. The verification results indicate a good agreement for temperature below 10 hPa and a relatively warm bias above. Since the warm bias between 0.4 and 10 hPa is similar for both MIPAS and HALOE, we conclude that the bias is due to the model or the meteorological analysis that drives GEM-BACH. Also, we should note that for this two-month period (1 August to 30 September 2003), the horizontal coverage of the two instruments is not the same. HALOE is limited to the band 40° S to 70° N for this time period.

The agreement in temperature below 10 hPa indicates that MIPAS temperatures are likely to be in agreement with radiosonde temperatures assimilated in GEM-BACH. Indeed, radiosondes measure temperature up to 30 hPa (in the tropics), but their effect on meteorological analyses can be observed up to about 10 hPa. Since radiosonde temperatures provide a strong constraint on meteorological analyses, they have a (significant) impact on GEM-BACH in the MR mode. Thus, the agreement with MIPAS temperatures below 10 hPa, as seen Figure 8, is an indication of agreement with radiosonde temperatures. Lastly, since we are specifically interested in temperatures over the polar region for heterogeneous chemistry, Figure S7 shows good agreement over the South Pole region (between 70° S and 90° S) during the month of August (important for the onset of ozone hole events).

7.2.2. Methane and Nitrous Oxide

Continuing our assessment of meteorology, evaluation of the distribution of long-lived species, in particular CH₄ and N₂O, can provide information about the quality of the wind fields. The evaluation of CH₄ against MIPAS and HALOE observations is presented in Figure 9. It shows that the model CH₄ is in good agreement with MIPAS observations across the entire stratosphere.

Also, remarkable agreement with MIPAS measurements of N₂O is observed (see Figure S8 in Supplementary Material). These results indicate that the wind fields produced by the meteorological analysis (in MR mode) are of good quality in the stratosphere. However, the comparison with HALOE CH₄ measurements shows non-negligible biases both in the lower and upper stratosphere, although the standard deviation is small, indicating that the spatial patterns in the modeled distribution of CH₄ remain close to the observations.

The fact that accurate winds can be obtained in the stratosphere deserves some attention. We recall that GEM-BACH is driven in MR mode. The meteorological analyses are affected by wind observations only from the troposphere. In the stratosphere, temperature-sensitive radiance observations are the main source of observations. Although satellite radiance observations often have an offset and may result in temperatures being inaccurate (even after the radiance bias correction), we can argue that the horizontal distribution of radiances is well captured, and consequently the horizontal temperature gradient is well represented. It is known that for the most part, on synoptic timescales, the stratosphere is in geostrophic balance (and, on large scales, in gradient-wind balance) as discussed in Section 2.2. Thus, we argue based on the thermal–wind relation Equation (3) that the vertical shear of the geostrophic wind in the stratosphere is also well captured since the horizontal gradient of temperatures is reasonably captured by the stratospheric meteorological analyses. The tropospheric winds are also well represented in tropospheric meteorological analysis. Thus, using these winds as the lower boundary condition in the thermal–wind Equation (3), we can deduce the 3D-distribution of the geostrophic wind in the whole stratosphere. This makes the point that the stratospheric winds are well represented. Lastly, since the wind field and tracers are closely related as a result of shear flow balanced by stirring and mixing, as discussed in Section 2.3, we can understand that the distribution of long-lived species is also well captured throughout the stratosphere with a coupled model driven in MR mode.

It is also useful to compare the verification by region (Tropics vs. South Pole) for the long-lived species CH₄ and N₂O. In an experiment where we gathered the data over the period of 17 August to 7 September 2003, we plotted the relative error (which is dimensionless, rather than the volume mixing ratio) for two regions, the tropics (from 30° S to 30° N) and over the South Pole (south of 60° S). The results are shown in Figure 10.

We note that above 3 hPa both modeled values for CH₄ and N₂O are overestimated in the tropics, and inversely, the modeled values are underestimated over the South pole region, whereas there is little bias below 3 hPa in the tropics. Since both species are chemical tracers, the effect presented is primarily due to transport, and is similar for each species. These results suggest that above 3 hPa the equator-to-south pole portion of the Brewer–Dobson circulation is too weak in the model, and thus compared with observations the modeled values are too small in the tropics and too large on the other extreme of the Brewer–Dobson circulation over the South Pole region.

7.2.3. Ozone and Water Vapor

Important meteorological-chemical interaction arises with gases such as O₃ and H₂O. Global averages of O₃ differences for the period of August-September 2003 are presented in Figure 11. We observe a significant model ozone deficit (i.e., underestimation) in the upper stratosphere with a maximum deficit of about 15% at 0.7 hPa. This can be explained, at least in part, by the model warm bias at these altitudes and the negative correlation between temperature and ozone, as explained in Section 2.1, in particular, making use of Equation (2). The model ozone deficit in the upper stratosphere may also be partly due the severe overabundance of model NO₂ in the upper stratosphere (see Figure S9 in Supplementary Material), since nitrogen dioxide, NO₂, catalytically destroys ozone. The model’s overestimation of NO₂ is also supported by the comparison of GEM-BACH against Fourier-Transform InfraRed (FTIR) spectrometer measurements at Eureka [109] (see also Section 7.3). A better agreement between GEM-BACH and HALOE in terms of O₃ is observed from 2 to 10 hPa, and with an overestimation below 10 hPa. This positive bias in GEM-BACH simulations is also seen when we compare the model against ozone sondes (displayed in Figure S10, Supplementary Material). The above considerations indicate that around 30 hPa GEM-BACH globally averaged O₃ concentrations are too high and that HALOE observations are more accurate than MIPAS ESA retrievals. However, over the South Pole region, Figure 12 shows that GEM-BACH concentrations match the MIPAS concentrations between 3 and 30 hPa, the region of ozone destruction.

Water vapor plays an important role in the ozone budget at the stratopause. It has an important radiative impact in the lower stratosphere/tropopause region. A comparison of GEM-BACH H₂O against MIPAS and HALOE profiles is presented in Figure 13. Good of agreement with MIPAS data throughout the stratosphere as well as a model overestimation near the model lid are evident.

We note that compared with HALOE data, there is a gradual increase in the H₂O bias with height, indicating a positive model bias which grows gradually with height across the stratosphere. Closer to the tropopause, we observe a significant model overestimation compared to both MIPAS and HALOE data. Since the retrieval error is important and the variability of H₂O is very large, a definitive conclusion would require a more detailed investigation, which we have not carried out here. In fact, since there is a sharp transition in H₂O at the tropopause and because the tropopause height varies considerably between the tropics and the polar regions, the global average is simply not a good comparison statistic.

7.3. GEM-BACH Evaluation against Ground-Based Total Column Measurements

An evaluation of stratospheric trace-gas column distribution in GEM-BACH was made against the Network for the Detection of Atmospheric Composition Change (NDACC) Bruker 125HR FTIR spectrometer data at Eureka (Eureka, Nunavut, Canada (80.05°N, 86.42°W)). The comparison was carried out during the International Polar Year (IPY), from 1 March 2007 to 28 February 2009 and the data are publicly available on the SPARC-IPY web site [110].

The FTIR instrument, retrieval methods and measurements during this time are described in depth in Batchelor et al. [111,112]. While an additional comparison between this FTIR and GEM-BACH and other models for the NOy budget has been published by Lindenmaier et al. [109], it is valuable to provide an example of how GEM-BACH compares with ground-based measurements as part of this discussion, especially for O₃, HCl and CH₄ which were not examined in [109].

Comparisons between the FTIR stratospheric partial columns measured at Eureka for six chemical species during 2007 are shown in Figure 14. Here the GEM-BACH model is refreshed with CMC meteorological 3D-Var-FGAT (First Guess at Appropriate Time) analyses every 12 h. The configuration of GEM-BACH is identical to that described in Section 3 except that the surface area density for the PSC’s has been reduced to provide a better agreement with ozone observations in polar regions.

As described in [111,112], the vertical resolution of the FTIR measurements is limited, and in some parts of the atmosphere the a priori contributes significantly to the column. To account for this and provide a clean comparison, the higher-resolution GEM profile ${\chi }_m$ is smoothed with the FTIR a priori ${\chi }_a$ and averaging kernel $\mathbf{A}$, to provide a smoothed profile ${\chi }_s$ as defined in [113],

$${\chi }_s={\chi }_a+\mathbf{A}({\chi }_m-{\chi }_a).$$

(19)

The column is determined from the smoothed profile, with the partial column bounded between 9.8 hPa (corresponding to approximately the bottom of the stratosphere) and the maximum altitude where data contributes more than the a priori (sensitivity > 0.5) to the FTIR retrieval. The argument to use Equation (19) can be understood from the fact that if we replace ${\chi }_m$ by the true concentration $\chi$, we get the standard expression for a retrieved profile using the averaging kernel $\mathbf{A}$. Then, with the ansatz that GEM-BACH provides an ideal high resolution profile, if we substitute it for the truth, we then get a smoothed or equivalent “retrieved” profile [113] that we can compare with the FTIR retrieval.

Observations at Eureka provide a usefully challenging test case for modeling, with the dynamic polar vortex allowing air inside, through the edge, and outside of the vortex to be sampled overhead. Figure 14 demonstrates how well GEM-BACH captures this dynamical variability. Between day 65 and 85 of the SPARC-IPY campaign, the FTIR sampled the air mass inside the polar vortex, as seen in the perturbed profile across all the gases. The 2006/2007 polar winter was characterized by a strong, cold polar vortex with significant amounts of ozone depletion [114,115]. This is well captured in the model, with ozone (panel (a)) tracking the FTIR columns extremely closely, both inside and outside the vortex. The day-to-day dynamic and seasonal variability is captured throughout the year across all gases, suggesting that both meteorological and radiation processes are being captured well. As described in [109] there are some consistent offsets seen between the model and data in the chlorine reservoir species HCl and ClONO₂, and in HNO₃ (though the latter matches well within the 2007 vortex, GEM-BACH is typically 10% lower than observed throughout the rest of the comparison period [113]). These differences are likely due to not including all of the chlorine sources (CFCs), as well as limitations in the PSC treatment, which is tuned to Antarctic conditions and does not include type 1b liquid PSC particles, which play a bigger role in the Arctic [113].

Methane (CH₄) is an important greenhouse gas, and several recent studies have focused on better understanding the contribution of the stratosphere to the methane column [116,117]. Figure 14 panel (f) shows a seasonal bias in the modeled stratospheric CH₄, with excellent agreement in the summer and a high bias in the spring and fall. A similar overestimation of CH₄ at high northern latitudes is observed and the cause remains a subject of investigation.

8. Summary and Conclusions

A fully coupled meteorology-chemistry model, called GEM-BACH, was developed by combining the stratospheric extension of the NWP model GEM with the BASCOE chemical transport model using a chemical interface that preserves time-splitting properties while being modular, allowing the system to run with or without chemistry. The increase in computational cost was minimal due to the semi-Lagrangian scheme where the upstream point and interpolation weights are computed only once for all species.

The project started with a preliminary stratospheric version of the CMC operational meteorological model GEM [57], where the non-orographic gravity wave drag [39,40] and a k-correlated radiative schemes [60] were added. The model then produced realistic lower stratosphere temperatures over the South Pole for the initiation of PSC chemistry, and with a prescribed Stratospheric Sulfate Aerosol (SSA) surface area densities gave rise to realistic simulation of the ozone hole without chemical data assimilation (but driven by meteorological analyses), a rare outcome for a pure simulation (Guy Brasseur, personal communication).

We compared chemical observations with GEM-BACH simulation in which with the meteorology component was replaced by meteorological analysis every 6 h. This mode of coupling with meteorological analyses was called meteorological refresh (MR) mode, and the accuracy of the simulated chemistry was compared with the standard CTM approach, where an offline chemistry model is driven by meteorological analyses linearly interpolated in time. The results show that the dynamical consistency provided by the coupled model driven with meteorological analysis refresh (MR), although with an offset created every 6 h by the meteorological analysis, is much more accurate than a linear interpolation of analyses used to drive an offline CTM. This conclusion was reached irrespectively of the type of meteorological analysis used, whether it came from a 3D-Var or a 4D-Var assimilation scheme.

We also add a discussion on the importance of having (relatively) high resolution meteorological analyses and analysis increments computed at the same resolution as in the coupled model. The GEM-BACH simulation of the ozone hole event of 2003, with the coupled model driven in MR mode, gave particularly good results compared with independent observations, in terms of values inside and outside the vortex, as well as the gradient along the vortex edge. The quality of the 4D-Var ozone assimilation performed with a CTM at much lower resolution did not approach the simulations with GEM-BACH, thus stressing the importance of resolution in obtaining accurate chemical fields.

The temperature in GEM-BACH, driven in MR mode with the standard meteorological observations showed fairly good agreement (in the lower stratosphere) with independent temperature measurements from MIPAS and HALOE. But large and similar biases against MIPAS and HALOE were observed in the middle and upper stratosphere, indicating a warm bias in either GEM or in the standard meteorological temperature data. On the other hand, the quality of the transport, evaluated by comparing the model output to observations of long-lived chemical species, showed good quality throughout the stratosphere, both in a global time-averaged mean and in their daily variability, compared with total column ground-based FTIR measurements. However, the model upper stratospheric ozone is underestimated by the BASCOE chemistry, which may be attributable to temperature overestimation or poorly modeled NOx at these altitudes. Rather good agreement was observed with chemically produced H₂O throughout the stratosphere.