## 1. Introduction

Most numerical models simulate momentum and energy fields reasonably well, but the question of how best to approximate continuity equations for trace constituents, which are fundamental in simulating climate (Chahine 1992) and effects of chemical constituents (e.g., Denning et al. 1995), has vexed modelers for decades. The underlying difference is that momentum and energy distributions remain reasonably continuous for atmospheric circulations so that the differential equations can be well approximated by finite-difference or spectral representations. In contrast, trace constituent distributions develop stronger horizontal and vertical gradients, resulting in substantial errors in finite-difference schemes and Gibbs oscillations in spectral models (Navarra et al. 1994). Consequently, a substantial variety of schemes has been developed to improve the accuracy of trace constituent transport computations. This paper describes experiments examining the accuracy of nine algorithms, using channel versions of the University of Wisconsin (UW) hybrid isentropic–sigma (*θ–σ*) model and the nominally identical sigma (*σ*) model.

Comparative tests of transport schemes have for the most part been either one- or two-dimensional. Under these limitations, transport schemes are often compared to a simple analytic solution under conditions in which a trace constituent distribution should retain its initial shape precisely throughout a simulation. In one-dimensional tests, initial shapes have included square blocks, triangles, and exponential “hills”; the most common two-dimensional distribution is probably a cone. Rood (1987) reviews such techniques and results of one-dimensional tests.

Relatively few efforts at comparing schemes have included the effects of the earth’s sphericity, the Coriolis force, and more importantly the three-dimensional effects of a vertically sheared atmosphere during the baroclinic amplification of a synoptic-scale wave. The experiments here compare the algorithms under these more difficult conditions. Although no exact analytic solution is available for comparison of results, three aspects of the solution are known. First, the algorithm should not develop negative trace constituent values. Second, the global integral of trace constituent should be exactly conserved; but this condition is forced in all the simulations and thus provides no means of comparing the accuracy of the schemes. Third, the experiments are all performed under adiabatic conditions, so that the trace constituent on any isentropic level should be conserved. This provides two means of comparing model accuracy:first, the maximum value of trace constituent on an isentropic surface, either explicit or implicit, should remain constant; and second, the area between any two contour values on an isentropic surface should remain constant.

The constraint of adiabatic conditions is justified in two ways. First, it provides two objective means of making intermodel comparisons under three-dimensional atmospheric motion. Second, for scales resolved by NWP and climate models, most trace constituent transport in the atmosphere occurs under conditions that are well approximated as isentropic. For example, Liu et al. (1995) determined that over the western equatorial Pacific Ocean warm pool only 15% of the clouds are precipitating at any one time.

The condition that a trace constituent maximum on an isentropic surface should remain constant also has implications beyond the present experiments. For example, precipitation occurs where water vapor achieves saturation; if the simulation of water vapor transport cannot accurately preserve its maximum value, condensation processes and precipitation will occur too soon or, more likely, too late or not at all (Johnson et al. 1993). Simulation of chemical processes requires the same degree of accuracy. Thus the condition for comparison of transport algorithms in these experiments is also an important requirement for all atmospheric models.

The condition that the area between any two trace constituent contour values should remain constant provides a more comprehensive basis of comparison than maximum values. These comparisons, suggested by an anonymous reviewer, provide some additional insights as to where particular models are more or less accurate, though by and large the conclusions are quite similar to those based on constituent maxima.

Tests of the conditions within a cyclic channel are with a standard *σ*-coordinate model with four different vertical resolutions, and the UW *θ–σ* model with two vertical resolutions. Because both models at all resolutions simulate standard meteorological quantities to a high degree of similarity through 48 h, they provide excellent test beds for insertion of various trace constituent transport algorithms. Further, the simulated results will be examined entirely within the *θ* domain of the *θ–σ* model; thus the experiments also provide comparisons of trace constituent transport algorithms as expressed in *θ* and *σ* vertical coordinates.

The results here are not intended to provide a comprehensive review of all possible transport algorithms; those included cover a range of sophistication. The nine schemes are standard second-order finite differencing (the “standard” scheme hereafter), the standard with three “borrow and fill” fixers, the “simple upstream” method (included only for completeness since its known dispersion has precluded its use in modeling efforts), the adaptation of Van Leer’s ideas (1974, 1977a, 1977b, 1979) by Allen et al. (1991), the adaptation of Van Leer by Lin et al. (1994), the “piecewise parabolic method” (PPM) as used by Carpenter et al. (1990), and the conservation of second-order moments (CSOM) scheme of Prather (1986). The simple upstream method, the Van Leer adaptations, and the PPM will be referred to collectively as the “conservative flux-integrated methods.”

The second section of this paper describes the models, how the inert trace constituent distributions are inserted into the initial atmosphere, and the transport algorithms. Section 3 presents the 48-h results, including the conservation of trace constituent maxima, conservation of contour areas, spatial distributions on 310 K, and vertical profiles. The fourth section reviews the results, implications, and limitations of the experiments.

## 2. The models, trace constituent distributions, and transport algorithms

### a. Description of the models and the initial atmosphere

The UW *θ–σ* channel model is hydrostatic and excludes orography. It is grid point in a scaled spherical domain extending from 10° to 70° meridionally and 45° zonally with east–west cyclic boundaries. Horizontal resolution is 1.5° in both directions. Walls at the north and south boundaries are impermeable and free slip. Vertically, the model has two domains: the sigma coordinate planetary boundary layer (PBL) comprising the lowest 100 mb; and an isentropic-coordinate free atmosphere extending from the top of the PBL to an upper boundary at 475 K. Two different models are used, labeled *θ–σ*(Δ*θ*), where Δ*θ* = 5 K or 3.33 K is the vertical resolution in the *θ* domain, except for the top layer where Δ*θ* = 120 K for both resolutions. The *σ* coordinate PBL has three evenly spaced predictive surfaces for both Δ*θ* resolutions. The predictive equations for mass and momentum utilize a second-order flux conservative algorithm, the “box” method (Kurihara and Holloway 1967). All simulations reported here are for 48 h, using a 3-min time step following Matsuno (1966). Additional details of model structure are found in Pierce et al. (1991) and Zapotocny et al. (1991, 1993).

Details of the *σ* model are provided in Johnson et al. (1993); it was adapted from the *θ–σ* model and is as close to identical as possible in all significant aspects with the exception of vertical coordinate. Results are reported for models of four different vertical resolutions denoted as *σ*(*n*) where *n* = 19, 27, 35, or 43 is the total number of vertical layers. For all *σ*(*n*), the PBL comprises three equally spaced layers of approximately 100 mb, duplicating as closely as possible the PBL of the *θ–σ* models.

The initial trace constituent distributions are inserted primarily above the PBL, and thus the relevant vertical resolution with respect to transport is in the free atmosphere. The *θ–σ*(5) model has an average of approximately 13 layers above the PBL (because isentropic grid volumes emerge and submerge, the number of layers is neither constant nor an integer), and the *θ–σ*(3.33) model has approximately 19 layers. The *σ*(19), *σ*(27), *σ*(35), and *σ*(43) models have 16, 24, 32, and 40 layers in the free atmosphere, respectively. Thus vertically the highest resolution *σ*(43) model has approximately twice as many layers as the highest resolution *θ–σ*(3.33) model in the domain of interest; the *θ–σ*(5) has the lowest resolution of all with approximately 13 layers.

The *θ–σ* initial atmosphere was generated with the analytic functions described in Pierce et al. (1991) using the parameters listed in Table 1 of Johnson et al. (1993). The *σ* initial atmosphere was interpolated from the *θ–σ* atmosphere as described in Johnson et al. After generation, all model atmospheres were dynamically initialized for 40 steps following Temperton (1976). All resolutions of both models produce 48-h predictions of mass and momentum, which are quite similar. Partial evidence of this statement is provided in Fig. 1: 500-mb geopotential heights for the *θ–σ*(5) model and the *σ*(43) model, in Figs 1a and 1b, respectively; and zonal velocities on 310 K for the same two models in Figs. 1c and 1d. (The models’ cyclic domain extends only 45° zonally; the area from 45° to 60° in Figs. 1, 4, and 5 is the same as from 0° to 15° and is included only for visual purposes.) The similarities in Fig. 1 are indicative of the results for all standard synoptic fields. (In the *σ* models, data on isentropic surfaces are obtained by interpolation with respect to *p*^{κ}, where *κ* = *R*_{d}/*c*_{p}.)

The 310-K layer remains entirely above the *σ* coordinate PBL throughout the 48-h integrations for both *θ–σ* models. Thus in comparisons of trace constituent distributions under isentropic conditions, differences between the *θ–σ* and *σ* models must be caused by differences, for each scheme, between the quasi-horizontal two-dimensional transport in the *θ* domain of the *θ–σ* models and the three-dimensional transport in the *σ* models. The results thus allow comparisons not only between the schemes but also between the attributes of *θ* versus *σ* coordinates. This question is addressed further in section 3.

### b. The initial trace constituent distributions

*q*

*H*

*V*

*q*represents an arbitrary trace constituent mixing ratio. For both distributions, the horizontal function

*H*is given by where

*λ, λ*

_{0}, and

*ϕ, ϕ*

_{0}(in degrees) specify meridional and zonal position, respectively. Here

*ϕ*

_{0}was chosen so that the distributions are centered in the channel zonally;

*λ*

_{0}is slightly equatorward of the channel center where the 310K isentropic surface is steeply sloped meridionally. Vertically the functions

*V*(

*θ*) and

*V*(

*p*) for the 310K and 660P distributions are, respectively, given by The pressure 660.6389 mb in (4) coincides exactly with an information surface of the

*σ*(43) model, and is also quite close to 310 K. The denominators 5 and 75 in (3) and (4), respectively, produce distributions that have about the same vertical extent with respect to pressure. Horizontal distributions of the initial trace constituent on the 310 K isentropic surface and north–south cross sections through the maxima are shown in Fig. 2.

Experiments were performed with two initial distributions in an attempt to provide “fair” comparisons between the *θ–σ* and *σ* models. Presumably the *θ–σ* models might enjoy an advantage (in some undefined sense) in simulating transport of the 310K distribution, which is centered vertically along a *θ–σ* model information surface; in turn, the *σ* models might have a similar advantage in simulating transport of the 660P distribution. The results in section 3 show that this concern was unfounded; for the most part, both the *θ–σ* and the *σ* models have more difficulty transporting the 660P distribution than the 310K.

### c. The transport algorithms

*t*can be symbolized as where

*q*is the constituent mixing ratio,

*u*is the zonal or meridional velocity component, Δ

*p*is the depth of the box, and g is the acceleration of gravity. Total trace constituent is automatically conserved, since any constituent transported out of one box is transported into its neighbor.

#### 1) The standard transport scheme

*i*and

*i*+ 1 is given by (5), where Comparable relations specify meridional and vertical transport.

This flux-conservative standard algorithm has many potential disadvantages, including the production of negative trace constituent values from “undershoots,” values larger than physically possible from “overshoots,” and the possibility of spurious oscillations as cited by Thuburn (1993). A further problem is that sharply “peaked” distributions tend to get dispersed. All the other algorithms utilized here were devised to avoid one or more of these problems.

#### 2) Borrow and fill techniques

Two common methods of controlling the development of negative trace constituents are borrowing and filling, used in three tests here, and the flux-corrected method of Zalesak (1979). In the methods here, a grid box with negative trace constituent is filled, ideally to a value of zero, and an equal amount is subtracted—borrowed—from one or more other boxes so that the model domain total constituent remains constant. Two such “Q-Fix” algorithms were used in experiments here, in which borrowing was from only the four nearest boxes (east, west, north, and south) within the same model layer.

After finding a grid box with negative trace constituent, the first Q-Fix algorithm (QF1) determines the neighbors with positive constituent and borrows an equal percentage from each. For example, if a box needs three units of trace constituent to be filled completely, and if the four neighbors have +10, +20, 0, and −1 units, respectively (+30 in the boxes with positive constituent), then 3 × (10/30) = 1 unit would be borrowed from the box with +10, 3 × (20/30) = +2 from the box with +20, and no constituent would be taken from the boxes with 0 or −1. Further borrowing was limited to one-third the total available. Thus, in the example above, if 15 units were needed with the four neighbors holding the same amounts as above, (⅓) × 30 = 10 would be borrowed (still in the ratio of one-third to two-thirds) and the box would be filled only to −15 + 10 = −5.

In the second Q-Fix method (QF2), borrowing was weighted by the amount available in neighboring boxes, as in QF1, and also by velocity components toward the boxes. QF2 extends the suggestion of Mahlman and Sinclair (1977), in their study of one-dimensional transport, that borrowing should first be from a neighboring box that is downstream. The idea can be understood by examining one way that negative values are created: if *q*_{i−1} = 0, *q*_{i} = 0, *q*_{i+1} > 0, and *u* > 0 (i.e., from box *i* toward *i* + 1), then the transport of (6) and (7) will produce a negative value for *q*_{i}. Under this condition, it seems reasonable that borrowing should be from box *i* + 1 and not from *i* − 1. As in QF1, borrowing was limited to one-third the total available, and if the amount available is insufficient, the box is not filled completely.

Neither QF1 nor QF2 completely removed all negative trace constituent. In order that results could be compared directly to the other algorithms, which by design eliminate all negatives, a “global” borrow and fill algorithm was implemented each time step after the“local” QF1 or QF2 algorithms. In the global fixer (following Rood 1987), the total negative (*N*) and positive (*P*) trace constituents (mass- and area-weighted) are computed. The negative trace constituent values are then replaced by zero, and each positive value *q*^{+} is replaced by (1 − *N*/*P*)*q*^{+}; this conserves the global integral of constituent.

In the interest of completeness, experiments were also performed with only the global fixer. For purposes of discussion below, the borrow and fill schemes (QF1 and QF2, with or without the global; and the global alone) will be referred to collectively as the “fixers.”

#### 3) Simple upstream

*q*in (5) is This simple upstream approach automatically eliminates negative values of trace constituent since if

*q*in the upstream box is zero, transport is also zero. This scheme is known to cause serious numerical diffusion.

#### 4) Two adaptations of Van leer

*q*

_{i}is a local extremum, a condition that prevents overshoots. Lin et al. (1994) compare this formulation with three others in one-dimensional transport of a rectangular wave. The experiments reported here use the choice of Lin et al. in their three-dimensional tests: the difference in

*q*at opposite sides of grid box

*i*is given by where

#### 5) Piecewise parabolic method

*δq*and

*a*

_{6}are parameters and

*x*is the normalized coordinate within a grid box: ½ ≥

*x*≥ −½. Similar to the slope in the two Van Leer algorithms, computation of

*δq*and

*a*

_{6}is undertaken by considering the values of 〈

*q*〉 in neighboring grid boxes. After

*δq*and

*a*

_{6}are computed, the function

*q*(

*x*) is examined to ascertain that it is positive and monotonic; if not, adjustments are made to the

*δq*and

*a*

_{6}values. Finally, the transported

*q*of (5) is given by ∫

*q*(

*x*)

*dx*where integration is over the

*u*Δ

*t*portion of the grid box, the cross-hatched area in Fig. 3b. Equation (11) has the attribute that the total constituent within a grid box is simply 〈

*q*〉Δ

*x*= 〈

*q*〉 because the integrals of the other two terms over the grid box are identically zero.

#### 6) Conservation of second-order moments

The last transport scheme is described by Prather (1986). The name Prather gives this scheme, conservation of second-order moments, is somewhat confusing but has been retained here because nothing else equally brief seems less so. This approach, like the conservative flux-integrated methods, assumes a polynomial expansion to describe the trace constituent distribution within each grid box. Ten orthogonal functions are used in the expansion, in which the functions comprise 10 terms: a constant; linear terms in *x, y,* and *η,* where *η* is the vertical coordinate *θ* or *σ*; quadratic terms *x*^{2}, *y*^{2}, and *η*^{2}; and cross terms *xy, xη,* and *yη.* The “moments” are the coefficients of these functions. In one sense, this method is an extension of the PPM algorithm in which the constant 〈*q*〉 and the functions *x* and ^{1}/_{12} − *x*^{2} are orthogonal over the domain of *x* from −^{1}/_{2} to ^{1}/_{2}. Further, as in PPM the constant term gives the total constituent within a box because the volume integral of the other functions is identically zero. However, the CSOM is substantially different in the computations and use of the moments. In the other schemes, the coefficients (i.e., the slope in the two Van Leer; and *δq* and *a*_{6} in PPM) are determined each time step from the *q*_{i} in a grid box and its neighbors; after the coefficients are used in the time step calculation of fluxes, they are discarded. As such, the only information retained in PPM for the next time step is the amount of trace constituent that crosses from one box to the next. The CSOM approach is more complicated: the spatial description of the trace constituent flux that crosses from one box to its neighbor, as embodied in the moments (the coefficients of the orthogonal functions), is also used in computing the spatial distribution for the next time step. A way of understanding the difference can be reached by examining Fig. 3b: in PPM, the only information transmitted from the box shown to its neighbor to the right is the *size* of the cross-hatched area, whereas in CSOM, the *shape* of the cross-hatched area (which varies in the directions perpendicular to the dimension shown in Fig. 3b, due to the *xy, xη,* and *yη* terms) is also transmitted. Thus the CSOM moments are used to compute both constituent averages and also distributions within the neighboring grid boxes. In contrast, the Van Leer and PPM coefficients are used to compute the neighboring constituent averages only.

Prather (1986) provides a computer listing for transport in one direction. The experiments here used the listing in three separate subroutines for zonal, meridional, and vertical transport, called sequentially. Experiments were undertaken to test whether the calling order of the three subroutines affected results: differences in the 48-h trace constituent distributions were very small, possibly no larger than those resulting from computer inaccuracies, and in no case large enough to produce visual changes in the distributions.

Unlike the conservative flux-integrated algorithms, the CSOM is not inherently positive definite. However, Prather (1986) discusses a “limit” procedure to prevent negative constituent from being transported from one grid box to another. When invoked in these simulations, no grid box developed negative trace constituent.

## 3. Results

### a. Negative trace constituent

As noted previously, the standard algorithm produces grid boxes with negative trace constituent. The ability of the QF1 and QF2 borrow and fill methods (without subsequent application of the global fixer) to eliminate negative trace constituent is indicated in Tables 1 and 2 for the 310K and 660P distributions, respectively. These tables show the total-domain negative trace constituent after 48 h as a percentage of the initial total positive. Since the standard algorithm conserves total trace constituent, the values in Tables 1 and 2 imply that the 48-h total positive constituent is greater than the initial.

As seen in Tables 1 and 2, QF1 and QF2 are approximately equivalent with respect to total-domain negative trace constituent; the small differences were not investigated. The CSOM scheme without Prather’s limit procedure also created negative values, though less than either QF1 or QF2. As expected, negative trace constituent did not occur in any experiment with the conservative flux-integrated schemes. It is also interesting to note that, for the standard case, both resolutions of the *θ–σ* model produce less negative constituent than any of the *σ* models for both distributions.

### b. Conservation of 310K tracer maximum

Under the adiabatic conditions of these experiments, the continuum equations require that the maximum value of trace constituent within each isentropic layer remain constant. Therefore deviations from the initial constituent maximum within explicit or implicit isentropic layers during integration is one measure of an algorithm’s accuracy. The tracer maxima were examined on 310 K because it is the most steeply sloped isentropic surface in the initial atmosphere, and thus provides the sharpest distinctions between models of different vertical resolution and coordinate.

Results of the 310K and 660P simulations are displayed in Tables 3 and 4, respectively, which provide the ratio of maximum trace constituent in the 310 K layer for the 48-h distributions to the maximum in the initial distribution. Under the adiabatic constraint, all values in Tables 3 and 4 should ideally be unity; the accuracy of trace constituent transport for a particular model, resolution, and/or scheme can be judged by the departure from this ideal. In these tables, the global fixer is employed after the two local schemes, QF1 and QF2, and Prather’s limit algorithm have been used. Thus all results in the tables, except for the standard method, are under conditions where negative trace constituent values have been eliminated.

Concentrating first on the two *θ–σ* columns in Tables 3 and 4, one aspect is immediately clear: the *θ–σ* model values are affected little by the change in vertical resolution. That is, regardless of fixer, transport scheme, or initial trace constituent distribution, the *θ–σ*(5) and *θ–σ*(3.33) values are quite close together. This result is easily understood: since flow in *θ* coordinates is two-dimensional under adiabatic conditions, differences between the *θ–σ*(5) and *θ–σ*(3.33) would be caused primarily by mass and momentum differences for the 310 K layers (not shown), which are smaller than the *θ–σ*(5) versus *σ*(43) differences seen in Fig. 1.

For the various *σ*(*n*) models, no consistent pattern in Tables 3 and 4 describes all the transport algorithms. For the conservative flux-integrated and CSOM schemes, there is a progression toward better conservation as vertical resolution increases, except for the 310K distribution for the PPM and CSOM. An increase in accuracy with increased model resolution is consistent with usual expectations. However, no comparable trend exists for the standard and fixer algorithms: as seen in Table 3 for the 310K distribution, the *σ*(19) model produces the value closest to unity and the value for *σ*(35) is the farthest, whereas in Table 4 for the 660P case, the *σ*(27) results are best and *σ*(19) and *σ*(43) are worst.

One objective in computing the data in Tables 3 and 4 is to compare the accuracy of conservation between the *θ–σ* and *σ* models. For the conservative flux-integrated and CSOM algorithms, the comparison is simple:for both initial trace constituent distributions, the *θ–σ* models have values closer to unity than the *σ* models. Only results from the CSOM *σ*(43) model are nearly as close to unity as those of the corresponding *θ–σ* models.

The comparison of *θ–σ* versus *σ* model results are less clear for the standard and three fixer schemes. If results for only the 660P initial distribution are considered, in Table 4, the *θ–σ* values are better: near unity for both *θ–σ* vertical resolutions but deviations as large as 10% for the *σ* models. However, this conclusion is contradicted by the results in Table 3 for the 310K initial distribution, where the *σ* values are better than those of the *θ–σ* models: the *σ*(19) model is best with less than 1% deviation from unity, followed by *σ*(27) and *σ*(43), with the *σ*(35) and two *θ–σ*(Δ*θ*) models the poorest and nearly the same at about 12% deviation.

Another stated purpose was to provide a means of judging which transport algorithm or fixer scheme is best for these particular three-dimensional experiments. The most consistent results, when considering both distributions and models, are obtained with the CSOM scheme. In particular, this scheme produces results distinctly the best in one respect: the values in Tables 3 and 4 are quite close to each other, for both the *θ–σ* and *σ* models; results for the other schemes depend more strongly on the initial distribution. However, Prather’s scheme as listed in his 1986 article increases computational requirements (using a scalar computer), for all models and resolutions, by about 20%. (As suggested by a reviewer, this disadvantage might be reduced if the scalar code were improved or a vectorized version used.) As a second conclusion, the standard second-order finite differencing plus either the QF1 or QF2 fixer is better than the conservative flux-integrated schemes, in spite of its dependence on initial distribution, for any *θ–σ* or *σ* model resolution.

Finally, these results have one completely consistent dependence on the initial distribution: every value in Table 4, for the 660P distribution, is less than the corresponding value in Table 3, for the 310K. The explanation is probably that the north–south extent of the 660P distribution (see Fig. 2) is less than in the 310K;or in other words, the 660P distribution is more sharply peaked. The dependence of maxima on initial distribution therefore suggests that all the schemes tend to disperse a more highly peaked distribution, even when the peaking is in only one dimension: the east–west extent of the two distributions is the same.

### c. Conservation of areas between trace constituent contours

As indicated in the introduction, not only should the maximum trace constituent value on any isentropic surface be conserved under the adiabatic constraint of these experiments, so also should the area encompassed between any two contours of constituent. Thus a comparison of these areas in the 48-h predictions to the corresponding areas in the initial atmosphere also provide a means of comparing the accuracy of different algorithms.

The data first extracted were the numbers of grid points between 0 and 1 tracer units, 1–2 units, and so on, except that values above the initial maximum of 10 were grouped with the 9–10 bin. [The maximum in the initial atmosphere is 10 tracer units, as indicated in Eq. (1).] The resultant 10 numbers for each of two initial constituent distributions, two *θ–σ* and four *σ* vertical models, and nine advection algorithms, proved to be somewhat overwhelming. Trial and error led to the conclusion that four larger bins were reasonably informative: between 0 and 1 units, 1–4, 4–7, and 7–maximum. The numbers of grid points in each bin were normalized against the corresponding values in the initial atmosphere, producing bias values, which in a perfect model would all have values identically 1. The resultant data for five of the algorithms and both initial atmospheres are shown in Table 5; data for the other algorithms are not included as they provide no information not obtained from the maximum value considerations of section 3b. The five algorithms included are the standard plus global fixer, Allen et al., Lin et al., PPM, and CSOM. As seen, it is quite difficult to summarize the data, but a couple of comparisons are of interest.

First, the standard plus global fixer “global only” algorithm does quite well in comparison to the more sophisticated Van Leer algorithms of Allen et al. and Lin et al. Values in the highest bin, 7–max, are clearly closer to unity and in agreement with the maximum value data of section 3b. For the other bins, the only place where the global only is not clearly better is in the 4–7 bin, for the 660P initial atmosphere; otherwise the data indicate no preference for the Van Leer algorithms, for either the *θ–σ* or *σ* models.

Second, the data on maximum trace constituent in section 3b suggested that the Allen et al. is preferable to the Lin et al. version of Van Leer. The more extensive data in Table 5 contradict this observation. It is true that the 7–max bin is preferable for the Allen et al., but for these two algorithms this bin is an extended version of examining the constituent maximum. For the 0–1 and 1–4 unit bins, the results provide no particular insight, but for the 4–7 bin, the Lin et al. algorithm appears somewhat preferable.

Third, the CSOM is once again the best of all the algorithms. Although a few of the bias scores for a particular bin, initial atmosphere, or algorithm at a specific resolution may contradict this statement, overall the CSOM clearly provides values nearest unity. This observation strengthens the conclusion formulated in section 3b.

### d. Horizontal distributions

Figures 4 and 5 show 48-h trace constituent distributions on 310 K. The results of standard second-order finite differencing are shown in Fig. 4 for three predictions with the 310K initial distribution: the *θ–σ*(5), *σ*(19), and *σ*(43). As with the 500-mb geopotential heights in Fig. 1, the *θ–σ*(3.33) results (not shown) are only slightly different than the *θ–σ*(5). The various *σ*(*n*) differences are somewhat more noticeable: as resolution increases from *σ*(19) to *σ*(43), the *σ* models transport progressively more tracer to the north side of the cyclonic circulation above the warm front. For example, the small eight-unit contour in *σ*(19) (Fig. 4b) becomes progressively larger, toward the south, north, and west in *σ*(43) (Fig. 4c). Distributions for the various borrow and fill fixers (not shown) are quite similar to those in Fig. 4: virtually no differences can be seen when only the global fixer is applied; and both local fixers, whether or not followed by the global fixer, produce a small intensification of the gradients along 40°.

Showing trace constituent distributions for both initial distributions, all transport methods, models, and vertical resolutions would require a great deal of space; little is lost by presenting only a few. Differences between the various transport schemes are seen in Fig. 5, of which five panels show results for the *σ*(43) 310K distribution, and one panel for the *θ–σ*(3.33). As anticipated, the simple upstream distribution (Fig. 5a) exhibits a severe decrease in the initial maximum. Judging the other algorithms in terms of the size of the highest value closed contour, the Lin et al. (1994) version (Fig. 5c) of Van Leer is somewhat less accurate than the Allen et al. (1991) (Fig 5b), and the PPM (Fig. 5d) of Carpenter et al. (1990), with its large closed 7 contour and small 8 contours, is somewhat more accurate than either. Prather’s (1986) CSOM (Fig 5e) is best since it retains a large 8 contour and isolated 9 contours.

The CSOM scheme from the *θ–σ*(3.33) model is shown (Fig. 5f) to make one comparison between the *θ–σ* and *σ* models: even though the *θ–σ*(3.33) model has only about half the number of layers of the *σ*(43) model in the free atmosphere, the 1.04 48-h constituent maximum of the *θ–σ*(3.33) model is slightly closer to unity than the 0.95 of the *σ*(43). This attests, at least for this case, to the ability of a model based on *θ* coordinates to simulate more accurate trace constituent transport than a significantly higher vertical resolution *σ* coordinate model.

### e. Vertical distributions

Thuburn (1993) cites a problem with second-order finite differencing: “Centered-difference schemes tend to produce spurious oscillations in the advected quantity upwind of sharp increases of gradient.” In the standard *σ* model experiments here, such oscillations were noticeable in vertical profiles at many locations for both initial atmospheres; an example for the 660P distribution is seen in Figs. 6a–d in the 48-h prediction at 43°N, 42°E. One instructive aspect of the profiles is that the ringing becomes progressively larger as vertical resolution increases. The reason for this behavior lies in the finite-difference approximations to the highly nonlinear vertical transport term, *σ̇**q*/∂*σ.* With baroclinic amplification of waves and the backing and veering of wind with height, higher-resolution *σ* models intensify both the constituent vertical gradient, ∂*q*/∂*σ,* and also the vertical motion *σ̇*

All four fixer routines (QF1 and QF2, with and without the global algorithm) remove the spurious vertical oscillations in the *σ* models at all resolutions; the results in Figs. 6e–h, for QF1 with global fixer, are typical of all the fixers. As in the horizontal distributions, the vertical profiles produced by the three fixers are virtually indistinguishable, for a given *σ*-model vertical resolution.

A purpose of Thuburn’s (1993) study was to demonstrate that a Van Leer approach, applied vertically in a spectral model, eliminates the oscillations. In fact, all four of the conservative flux-integrated schemes and the CSOM do so nicely in all resolutions of the *σ* model. Results are shown for the PPM in Figs. 6i–l; the other flux-integrated and CSOM schemes achieve similar results, though the maximum values in the profiles at a given location are usually different.

Profiles for the standard versions of the *θ–σ* models also show negative values; but since the simulations are under isentropic conditions, they are caused only by horizontal finite difference computations and thus, as evident from the examples in Figs. 6m,n, the *θ–σ*(5) and *θ–σ*(3.33) negatives are about the same size. Ringing in the *σ* model vertical profiles cannot, and does not, occur in the *θ–σ* model from vertical transport. All the fixers applied in the *θ–σ* models eliminate the negatives in much the same way; Figs. 6o,p show the results of fixer QF1.

Finally, increases in *σ* model vertical resolution reduce vertical diffusion in all four conservative flux-integrated and CSOM schemes. This statement is verified by data on the vertical width of constituent profiles, defined as the difference between the two pressures (below and above) at which the trace constituent is *e*^{−1} times its maximum in a profile. The vertical widths for all *σ*(*n*) models of the 48-h 310K and 660P distributions, averaged over 77 profiles, are shown in Tables 6 and 7, respectively. (The profiles are all located within the two-units contour of Fig. 4c.) As seen, for all five algorithms and for both initial distributions, the average vertical width of the distribution in the 48-h simulations becomes progressively smaller as *σ* model vertical resolution increases. The CSOM algorithm also maintains the narrowest vertical extent for nearly all resolutions and algorithms examined, when comparing *σ* models of the same resolution.

As pointed out by a reviewer, the ringing caused by second-order finite differencing could also be controlled by selective vertical filters. Indeed, the fixers are a form of filter that reduce the formal accuracy of the computations to first order. The CSOM limit procedure also reduces the formal accuracy of the scheme; but since the amounts of trace constituent to be eliminated (see Tables 1 and 2) were at least four orders of magnitude smaller than that generated by standard second-order finite differencing, this is a technical point and probably of little practical significance.

## 4. Discussion

The experiments presented in this paper examined the ability of nine transport algorithms to conserve the initial maximum on an isentropic surface within a baroclinically amplifying synoptic-scale wave. The results suggest several conclusions. First, standard second-order finite differencing creates ringing in the *σ* model vertical profiles from vertical transport, which cannot occur under adiabatic conditions in the *θ–σ* model. Further, increases in *σ* model vertical resolution cause intensified ringing. The question of how best to formulate a “borrow and fill” fixer to eliminate these distortions is more difficult for a *σ* model than for the *θ* domain of the *θ–σ* model. The apparent similarity in application of the QF1 and QF2 algorithms to the *σ* and *θ–σ* models masks a fundamental difference. The similarity is that both schemes borrow from the nearest neighbors on the same model layer, regardless of *σ* or *θ–σ.* However, the transport creating the negatives is very different. In the *θ* domain of the *θ–σ* model, negative values are caused entirely by horizontal transport, whereas the negative values in the *σ* models are caused by both horizontal and vertical transport. Thus a degree of consistency exists between the causes and subsequent elimination of negatives in fixers applied in the *θ* domain of the *θ–σ* model: the negatives are caused by excess removal of constituent in the calculation of horizontal transport between horizontal neighbors, and thus the nearest horizontal neighbors are the logical place to borrow. There is no easily implemented means of achieving a similar consistency for borrowing in a *σ* model.

A second conclusion is that the accuracy of a particular transport scheme must be tested three-dimensionally. Lin et al. (1994) provide an argument that their application of Van Leer’s principles is “less diffusive and therefore more accurate” than that of Allen et al. (1991), the other Van Leer scheme examined here. As seen in Tables 3 and 4, the Lin et al. scheme produced slightly larger maxima on the 310 K level when applied within the two-dimensional quasi-horizontal transport of the *θ–σ* models; but when applied in the *σ* models, the Allen et al. approach produced better maxima in these experiments—distinctly better for the 660P initial distribution. On the other hand, Lin et al. appears slightly preferable in the bias scores of the conservation between constituent contours data of section 3c. Finally, both these Van Leer extensions were formulated as the result of perceived weaknesses in second-order finite differencing; these weaknesses are not apparent in Table 5.

Third, the *σ* models clearly produce smaller maxima on 310 K than occurs within the *θ–σ* models, for any of the four conservative flux-integrated algorithms. The difference in maxima between the *σ* and *θ–σ* models is also significant for Prather’s CSOM scheme; only the highest resolution *σ*(43) model, with approximately twice as many vertical layers as the *θ–σ*(3.33), produces maxima nearly as accurate.

Fourth, the experiments here indicate that, despite its relative complexity and computational inefficiency (about 20% more computer time compared to the fixers, using Prather’s 1986 scalar listing), the CSOM algorithm is worthy of wider experimentation. Compared to the other algorithms, differences in trace constituent maxima (Tables 3 and 4) between the 310K and 660P distributions are smallest for the CSOM scheme. Differences between the 310K and 660P distributions are also smaller for the CSOM versus other algorithms for most of the bins presented in Table 5, particularly the 1–4 and 7–max bins. The CSOM algorithm also maintains the narrowest vertical extent of each distribution (Tables 6 and 7).

Fifth, in terms of negative constituent and conserved maxima on 310 K, these experiments are not conclusive with respect to the question of whether standard second-order finite differencing in the *θ–σ* or *σ* models, with or without fixers, produces the most accurate 48-h predictions. The *θ–σ* models produce only about two-thirds the negative trace constituent of the *σ* models, but the QF1 and QF2 fixers applied in the highest resolution *σ*(35) and *σ*(43) models produce the least negative (see Tables 1 and 2). Further, the *σ* models conserve trace constituent maxima slightly better for the 310K distribution (Table 3), but the *θ–σ* models are distinctly better for the 660P distribution (Table 4). In addition, there is no consistent trend of accuracy with vertical resolution in the *σ* models. No simple explanation seems to exist for these contradictory results. However, the question as to whether a *θ* or *σ* model produces the most accurate constituent maxima has been addressed by Zapotocny et al. (1997), in 10-day integrations using a global version of the UW *θ–σ* model and assimilated data for initial conditions. In these longer simulations, the trace constituent develops sharper vertical gradients that are readily simulated by the *θ* domain of the UW *θ–σ* model but with considerably less accuracy by the comparable UW *σ* model and the Community Climate Model 2 of the National Center for Atmospheric Research.

The experiments here do not test the question of whether the *θ* coordinate is better when diabatic processes are included. Comparisons of comparable *θ–σ* and *σ* models with diabatic processes have been completed (Johnson et al. 1993; Zapotocny et al. 1997), but not for the other algorithms tested here. An important conclusion from these results is that standardized tests of numerical schemes, as suggested by Pielke et al. (1995), should include the choice of vertical coordinate and be carried out three-dimensionally.

Semi-Lagrangian transport (SLT) has recently become a subject of interest. Part II of this paper (Reames and Zapotocny 1999) examines a variety of SLT algorithms using the same *θ–σ* and *σ* channel models and comparison techniques employed here.

## Acknowledgments

Support of the Department of Energy under Grant DE-FG02-92ER61439 and of NASA under Grants NAG5-1330 and NAG5-4398 is gratefully acknowledged. Thanks are due to Prof. Donald R. Johnson for careful reading of the manuscript and his many valuable comments, Mr. Todd K. Schaack for his gracious input, Judy Mohr for technical typing assistance, and Krista Ommodt for preparing the diagrams. We also thank our anonymous reviewers for suggesting the methodology of comparing areas between trace constituent contours as well as other helpful input.

## REFERENCES

Allen, D. J., A. R. Douglass, R. B. Rood, and P. D. Guthrie, 1991: Application of a monotonic upstream-biased transport scheme to three-dimensional constituent transport calculations.

*Mon. Wea. Rev.,***119,**2456–2464.Carpenter, R. L., Jr., K. K. Droegemeier, P. R. Woodward, and C. E. Hane, 1990: Application of the piecewise parabolic method (PPM) to meteorological modeling.

*Mon. Wea. Rev.,***118,**586–612.Chahine, M. T., 1992: The hydrological cycle and its influence on climate.

*Nature,***359,**373–380.Denning, A. C., I. Y. Fung, and D. Randall, 1995: Latitudinal gradient of atmospheric CO2 due to seasonal exchange with land biota.

*Nature,***376,**240–243.Johnson, D. R., T. H. Zapotocny, F. M. Reames, B. J. Wolf, and R. B. Pierce, 1993: A comparison of simulated precipitation by hybrid isentropic–sigma and sigma models.

*Mon. Wea. Rev.,***121,**2088–2114.Kurihara, V., and J. L. Holloway, 1967: Numerical integration of a nine level global primitive equation model formulated by the box method.

*Mon. Wea. Rev.,***95,**509–529.Lin, S.-J., W. C. Chao, Y. C. Sud, and G. K. Walker, 1994: A class of the Van Leer-type transport schemes and its application to the moisture transport in a general circulation model.

*Mon. Wea. Rev.,***122,**1575–1593.Liu, G., J. A. Curry, and R.-S. Sheu, 1995: Classification of clouds over the western equatorial Pacific Ocean using combined infrared and microwave satellite data.

*J. Geophys. Res.,***100,**13 811–13 826.Mahlman, J. D., and R. W. Sinclair, 1977: Tests of various numerical algorithms applied to a simpler trace constituent air transport problem.

*Fate of Pollutants in the Air and Water Environments,*I. H. Suffet, Ed., Vol. 8, John Wiley, 223–252.Matsuno, T., 1966: Numerical integrations of the primitive equations by a simulated backward difference method.

*J. Meteor. Soc. Japan,***44,**768–784.Pielke, R. A., and Coauthors, 1995: Standardized test to evaluate numerical weather prediction algorithms.

*Bull. Amer. Meteor. Soc.,***76,**46–48.Pierce, R. B., D. R. Johnson, F. M. Reames, T. H. Zapotocny, and B. J. Wolf, 1991: Numerical investigations with a hybrid isentropic–sigma model. Part I: Normal-mode characteristics.

*J. Atmos. Sci.,***48,**2005–2024.Prather, M. J., 1986: Numerical advection by conservation of second-order moments.

*J. Geophys. Res.,***91,**6671–6681.Reames, F. M., and T. H. Zapotocny, 1999: Inert trace constituent transport in sigma and hybrid isentropic–sigma models. Part II:Twelve semi-Lagrangian algorithms.

*Mon. Wea. Rev.,***127,**188–200.Rood, R. B., 1987: Numerical advection algorithms and their role in atmospheric transport and chemistry models.

*Rev. Geophys.,***25,**71–100.Temperton, C., 1976: Dynamic initialization for barotropic and multi-level models.

*Quart. J. Roy. Meteor. Soc.,***102,**297–311.Thuburn, J., 1993: Use of a flux-limited scheme for vertical advection in a GCM.

*Quart. J. Roy. Meteor. Soc.,***119,**469–487.Van Leer, B., 1974: Towards the ultimate conservative difference scheme. Part II: Monotonicity and conservation combined in a second-order scheme.

*J. Comput. Phys.,***14,**361–370.——, 1977a: Toward the ultimate conservative difference scheme. Part III: Upstream-centered finite-difference schemes for ideal incompressible flow.

*J. Comput. Phys.,***23,**263–275.——, 1977b: Toward the ultimate conservative difference scheme. Part IV: A new approach to numerical convection.

*J. Comput. Phys.,***23,**276–299.——, 1979: Toward the ultimate conservative difference scheme. Part V: A second order sequel to Godunov’s method.

*J. Comput. Phys.,***32,**101–136.Zalesak, S. T., 1979: Fully multidimensional flux-corrected transport for fluids.

*J. Comput. Phys.,***31,**335–362.Zapotocny, T. H., D. R. Johnson, F. M. Reames, R. B. Pierce, and B. J. Wolf, 1991: Numerical investigations with a hybrid isentropic–sigma model. Part II: The inclusion of moist processes.

*J. Atmos. Sci.,***48,**2025–2043.——, ——, and ——, 1993: A comparison of regional isentropic–sigma and sigma model simulations of the January 1979 Chicago blizzard.

*Mon. Wea. Rev.,***121,**2115–2135.——, A. J. Lenzen, D. R. Johnson, T. K. Schaack, and F. M. Reames, 1997: A comparison of inert trace constituent transport between the University of Wisconsin isentropic–sigma model and the NCAR Community Climate Model.

*Mon. Wea. Rev.,***125,**120–142.

Total negative trace constituent at 48-h for the 310-K distribution: percentage of 00-h trace constituent.

Total negative trace constituent at 48-h for the 660P distribution: percentage of 00-h trace constituent.

Conservation of maximum trace constituent on the 310-K surface for the 310K distribution: *Q*_{max}(48 h)/*Q*_{max}(0 h).

Conservation of maximum trace constituent on the 310-K surface for the 660P distribution: *Q*_{max}(48 h)/*Q*_{max}(0 h).

Conservation of areas between contours: ratios of 48-h predictions to values in the initial distributions.

Vertical width (mb) of 48-h 310K distribution averaged over 77 profiles.

Vertical width (mb) of 48-h 660P distribution averaged over 77 profiles.