Gridsize induced error in the discretization of exchange
processes at the tropopause

 

M. FIEBIG-WITTMAACK
Departamento de Matemáticas, Universidad de La Serena
Centro de Estudios Avanzados en Zonas Áridas, Casilla 599, La Serena, Chile.
E-mail: mefiebig@userena.cl

 

Received January 26, 2005, accepted April 26, 2005

 

RESUMEN

Estudiamos el error introducido por el método de las diferencias finitas en la discretización de un modelo global simplificado 2-D de transporte de gases traza, para casos en que los coeficientes de difusión, que relacionan el flujo con el gradiente de la razón de mezcla, tienen discontinuidades de salto en la tropopausa. Analizamos el método convencional de celdas tanto para el caso de un flujo ascendente típico, como también para el caso de un flujo descendente típico con reacciones químicas, comparando las aproximaciones de las soluciones correspondientes para diferentes tamaños de paso de discretización. Para el flujo descendente típico resulta que si la rejilla no es suficientemente fina se pueden generar grandes errores; estos se propagan fundamentalmente en la troposfera. En cambio, el flujo típicamente ascendente resulta ser relativamente insensible al tamaño de paso de la discretización.

 

ABSTRACT

We study the accuracy of the finite differences discretization scheme for a 2-D simplified model of global tracer transport, in the case that the diffusion coefficients relating flux to the gradient of the mixing ratio have discontinuity jumps at the tropopause. We analyze the conventional box method for a typical downward flow with chemical reaction and for a typical upward flow, comparing the approximations of the solutions, for different discretization gridsizes. It turns out that the jumps may introduce remarkable errors in the discrete solutions, in the case of a typical downward flow; these errors propagate mainly into the troposphere. A noticeable improvement is achieved by reducing the gridsize. However, a typical upward flow is rather insensitive to the chosen gridsizes.

Keywords: Global tracer transport model, finite differences discretization scheme, tropopause discontinuity jumps, discretization errors.

 

1. Introduction

The earth's atmosphere is divided into several layers. Of particular interest for the global transport of tracers are the troposphere (from ground up to about 8 km height near the poles, about 18 km height near the equator) and the stratosphere (the next higher level). There results a rather sharp boundary between the two layers, the so-called tropopause, below of which vertical exchange coefficients are large, whereas above it they are small (Birner et al., 2002; 2003; Schneider, 2004; Wirth, 2004).

Therefore, in 1-D and 2-D models of global tracer transport, the eddy diffusion coefficient tensor, relating flux to the gradient of the mixing ratio of some tracer has a discontinuity jump at the tropopause and also contains non-diagonal matrix elements in the lower stratosphere (Gidel et al., 1983; Brühl and Crutzen, 1988). The question arises whether in numerical models a significant error is introduced by the numerical schemes, due to the tropopause jump. This question seems to be especially relevant if we consider the increasing interest in recent years of the atmospheric science community to understand stratosphere-troposphere exchange processes and also, to understand the impact of aircraft emissions on atmospheric chemistry. Both issues were analyzed in many model studies (Ko and Douglass, 1993; Brühl et al., 1998; Groâ et al., 1998; Rondanelli etal., 2002; Gallardo et al., 2004); in almost all of them the effect of the jump discontinuity of the coefficients, on the accuracy of the numerical scheme, is not considered.

In this paper we show the error behavior of the conventional box method, for the 2-D case (Gidel et al., 1983; Brühl, 1987), with simplified numerical examples.

The basic equation for the time evolution of the concentration X of a particular tracer species in the atmosphere is given by

.....................................................(1)

where F is the flux density vector of this quantity and Q is its net chemical production per volume unit. With u = X/M_air denoting the mixing ratio of this tracer, M_air being the air concentration,

...........................................(2)

consists of an advective term with being the velocity vector of an averaged motion, and of a second term describing eddy diffusion with IK being the symmetric positive definite tensor of eddy diffusion coefficients.

Since our emphasis is on the effects of the jump in the diffusion coefficients with the conventional box method, we adopt a very simple model which only mimics diffusion and a chemical reaction, thereby neglecting advection and the effects of spherical geometry. Also spatial variations in the coefficients are not taken into account, besides of the variations resulting from the preferentially isentropic diffusion in the stratosphere which is closely connected to the shape of the tropopause. This does not mean that the neglected effects are not significant, but we want to clarify the numerical effects of the jumps without mixing in many different phenomena.

We restrict the chemistry term to a simple first order reaction:

.............................................................(3)

with constant reaction rate . Furthermore we treat the air density as a constant such that equation (1) reduces to

....................................(4)

In addition, suitable boundary conditions have to be stated.

 

2. The simplified 2-D case

We consider the continuous model for a simplified cartesian case without advection and with constant air density. The corresponding differential equation for the mixing ratio u of a tracer species is

.....................(5)

with

..................................(6)

and initial condition

.............................(7)

K, L and M are the space dependent exchange coefficients, i.e. the "eddy diffusion coefficients", and satisfy the condition

..............................................................(8)

for the sake of positive definiteness of IK with

The chemical reaction rate is considered constant. We suppose that at the lateral boundary the flux is purely vertical, i.e. the normal flux F_y is zero; but since we further assume that the off-diagonal diffusion coefficient L is zero at the boundaries we have

..................(9)

at the lateral boundary. Therefore the boundary conditions are given by (9) and either in the case

..................(10)

of tracers characterized by an "upward flow" (e.g. the CFC-family), or in the case of tracers

..................(11)

characterized by a "downward flow" (e.g. the ClO_x family). The matching conditions at the jump line are well known (Marchuk and Skiba, 1976; 1992), i.e. at the tropopause these conditions are now given by

................................................................(12)

and

.........(13)

where is the slope of the tropopause and the notations and denote the limit values of w taken from the stratospheric and the tropospheric side, respectively. Condition (12) is a continuity condition for the mixing ratio u at the tropopause and condition (13) is a continuity condition for the transverse flux component. We further suppose that the altitude z_TP of the tropopause is a piecewise continuously differentiable function of y, i.e. z_TP = (y) for 0< y < Y.

 

3. The conventional box discretization

Several investigations in air chemistry (McRae et al., 1982; Gidel et al., 1983; Brühl, 1987) solve equation (1) numerically by applying a splitting method (Marchuk, 1975) where the transport step uses Euler's forward method as time discretization and for concentrations and fluxes centered difference methods are used as spatial discretization (Fiebig-Wittmaack and Börsch-Supan, 1994). This conventional method uses discretizations based on cells, or boxes, defined by the inequalities

(j-i) Δy < y < j Δy and (i-1) Δz < z < iΔz, i=1,...,n; j=1,...,m,

with Δy = Y/m and Δz = Z/n. The grid G is given by

G= {(y_j,z_i)/y_j = (j-1/2) Δy, z_i=(i-1/2) Δz, i=1,...,n; j=1,...,m}

and the discrete values u_j,i of the mixing ratio are given at the midpoint of the cell. In order to improve the discrete boundary conditions in our simplified case, we introduce additional fictive outer points. We denote by F_y^j,i the discrete value of F_y calculated at the midpoint of the right boundary of j, i- th cell and by F_z^j,i the discrete value of F_z calculated at the midpoint of the upper boundary of the j, i-th cell. Therefore, we obtain the following discrete equation for the j, i-th cell:

..............(14)

The values at the fictive outer points must be calculated at every time step from the discretized boundary conditions. In order to ensure D-stability the inequality must hold; otherwise small errors

..............(15)

would propagate exponentially in time, the faster the smaller the time step is.

 

4. Numerical example

We use simplified data taken from the distributions employed in the 2-D box model of the MPI Mainz (Brühl, 1987; Fiebig-Wittmaack and Börsch-Supan,1994). Because of symmetry with respect to the equator we consider only one half sphere with y = 0 corresponding to the equator and y = Y = 9990 km corresponding to the north pole; z = 0 and z = Z correspond to altitudes of 2 km altitude above ground (just above the planetary boundary layer) and 23.6 km altitude (near the ozone layer), respectively. We introduce the altitude above sea level = z + 2. In agreement with the MPI model our basic gridsizes are Δz = 1.8 km, Δy = 1110 km (i.e. for the latitude φ we take Δφ = 10° ), and the time step is Δt = 2 hours.The eddy diffusion coefficients are given by

K=6x10³ Km²/h

with

Here may be interpreted as the slope of the main diffusion direction; thus this deviates from the horizontal in the stratosphere between 25° and 75° latitude. Note that a = 0 near the upper boundary has been introduced artificially in order to simplify the upper boundary conditions.

The tropopause is approximated by the following description

 

The shape of the tropopause and the eddy diffusion coefficients L and M are plotted in Figure 1 and Figure 2 respectively.

We concentrate in two cases, namely an upward flow with = 0 and the boundary conditions (9) and (10) with k_z = 5xlO^-4km/h (modeling CFC flow) and a downward flow with = 4x10^-4h^-1 and boundary conditions (9) and (11) with the "absorption" coefficient atthe lower boundary k_b=-3.5xl0^-3 km/h (modeling ClO_x flow). In both cases the prescribed boundary flux density is normalized to ± 1.

Stationarity is achieved by applying (14) over a sufficiently large time interval, until the computed values do no longer change significantly in time (nearly 19,000 iterations were needed in the " CFC case" and nearly 30,000 iterations were needed in the "ClO_x case"). We shortly describe how the discrete solutions look like: In the "CFC case", along vertical lines, the mixing ratio decays upwards in a nearly piecewise linear fashion, from values near 9 units at the lower boundary to values near 2 units at the upper boundary, with rather small differential quotients in the troposphere (0.05 units per km near the equator to 0.04 units near the pole), and larger ones in the stratosphere (from 0.8 to 0.5 units per km); along horizontal levels one sees a maximum of the mixing ratio at the tropopause. In the "ClO_x case" there is an exponential decay downwards of the tracer's mixing ratio by factors between 0.5 per km (near the equator) and 0.6 per km (near the pole). In the troposphere the decay is much smaller (from 0.85 to 0.9 per km); along horizontal levels one sees a minimum of the mixing ratio at the tropopause.

Since the exact solutions are not known, we improve the numerical approximations by reducing the step sizes to Δy/3 and Δz/3 first and to Δy/9 and Δz/9 then, and estimate errors by comparison to the next more accurate case. Differences in percentage from the discrete solutions for the downward flow are shown in Figure 3 and Figure 4, and for the upward flow in Figure 5.

 

In observing the errors shown in Figures 3, 4 and 5, it should be kept in mind that these are global errors as propagated from their source points to the point considered in the figure and do not simply reflect the local errors at the points of the grid. In the 1-D case, by evaluating Green's function and relating it to the solution, in the "ClO_x case", one can show also theoretically (see Appendix), that locally generated relative errors are propagated essentially downwards whereas the influence upwards decreases exponentially. In the "CFC case" the relative influence is rather constant downwards and weakly decreasing upwards. We observe a similar behavior in the 2-D case.

A consequence of this behavior is that, in spite of the fact that relative local errors are much larger in the stratosphere than in the troposphere because of larger variability of the solution there, the maximum of relative global errors lies in the troposphere where all errors propagated into.

 

5. Conclusions

The examples treated with our simplified model suggest that an upward flow without chemical reactions as typical for CFC, is rather insensitive to the large gridsize of the conventional treatment of the tropopause region (Fig. 5).

For a downward flow with reaction, as typical for ClO_x, large relative errors will occur already in the lower stratosphere (Fig. 3), since the equilibrium between diffusion and chemical destruction leads to a nearly exponential distribution of the mixing ratio in which the values in vertically neighboring cells differ by a factor between 2 and 4, in our example. If this is to be avoided, shortening of the meshsize is necessary (Fig. 4).

We conclude that for some tracers the jump of the coefficients at the tropopause can introduce remarkable errors in the discretizations, which propagate mainly into the troposphere.

In order to study the exchange processes at the tropopause, or, the impact of aircraft emissions, a careful numerical treatment of the jump effects near the tropopause must be achieved in the models; either shortening of the meshsize, or, a completely consistent treatment as in the Ritz-Galerkin finite element solution, eventually with triangular elements, may be the better way.

Finally, we conclude that even if we agree that important uncertainties of air chemistry models stem from the performance of the underlying dynamical model (Austin et al., 2003), it must also be kept in mind, that the numerical treatment itself can introduce significant errors in the models output if the gridsize is not chosen adequately.

Further work should treat extensions of our simple model by introducing advection with a realistic meridional circulation, including spherical geometry and variable air density. Also the spatial-seasonal variation of eddy diffusion has to be taken into account in order to achieve a realistic picture.

 

Acknowledgments

We thank Professor Dr. Wolfgang Boersch-Supan who suggested the subject and contribute with helpful discussions to its progress. Comments on the original manuscript by two anonymous reviewers are greatly appreciated. This work was supported by CONICYT-Chile under Grant FONDECYT 1030809.

 

Appendix: The 1-D case

The only dimension considered here is the altitude z. The vertical exchange coefficient M is called K in this case. The subscript S refers to the stratosphere, whereas T refers to the troposphere. The lower boundary of the domain is z = 0 or = 2 km respectively, the tropopause is assumed to be at at or , and the upper boundary of the domain is at z = Z = 21.6 km or z = 23.6 km. The values used for the exchange coefficientes are K_S = 1.5 x 10^-3 km²/h and K_T = 3x10^-2km²/h.

We further introduce the thicknesses of the stratosphere and of the troposphere: H_S = Z - z_Tp = 12.6 km and H_T = z_Tp = 9 km. The chemical reaction constant is = 4 x 10^-4/h in the ClO_x case (as with 2-D ), and = 0 in the CFC case. The upper and lower boundary conditions are the same as in the 2-D case. From the differential equation there result the exponential coefficients for the solution: . The boundary value problem can be solved exactly.

For the mixing ratio u we obtain in the "ClO_x case" the following solution:

The constants are given by A_S = 1.938, B_S= -1.231, A_T = 0.707, B_T = - 4.42 X 10^-4. At the upper boundary, the B_S-term can be neglected and we obtain u(Z) = 1.291 x 10³. Going downward we observe a fast decline of the solution corresponding to the A_S-term until near the tropopause where the B_S-term comes into play. At the tropopause we have u(z_Tp) = 0.707, and then a slow decline until we reach u(0) = 0.248 at the lower boundary.

We now proceed to the Green's influence function G(z ; ) from a unit source at . If the source is in the troposphere (i.e. < 0) we obtain

with

Due to the dependence on a further simplification is not possible. Note, however that 1 > exp > 0.353 and 1<exp< 2.83.

We observe that at the upper boundary the two terms are equal. Going downwards the B_S-term dominates and grows down to the tropopause. Also in the troposphere above the source point, the BTo-term dominates the A_To-term such that we have further growth there. Below the source point, however, the A_Tu-term is dominating causing a decrease of the function G downwards. This shows that the influence of an error in the troposphere relative to the solution remains about constant in the downwards direction whereas it declines in the upwards direction.

If the source is in the stratosphere ( i.e. > 0 ) we obtain

with

Again, a further simplification is not possible, but note that

We have a similar behavior as in the case of the source point in the troposphere, but between the source point and the tropopause the function G is decreasing downwards. The statement about the relative error remains unchanged.

In the special case = 0 we can simplify:

with

This function declines from the tropopause height z_Tp to both sides exponentially. Relative to the solution of the original problem this means that the influence of an error at the tropopause remains constant down to the lower boundary whereas it declines very fast in the upward direction.

In the "CFC case", because of missing chemical reaction, the solution is piecewise linear. We have

From the simple differential equation, and the boundary and matching conditions we obtain

and hence,

The Green's function for this problem is treated with a similar ansatz, however splitted at the source point:

If >0then

with

hence

Relative to the above solution one sees that below the source point the influence of an error is reduced, the more the larger , whereas above the source point it remains constant.

If <0 then

with

hence

Relative errors behave analogously to the case > 0, but the dependence on is weaker. Eventualy, for z = 0, we obtain

with

hence

We see that the relative error remains nearly constant over the whole domain.

 

References

Austin, J., D. Shindell, S. R. Beagley, C.Brühl, M. Dameris, E. Manzini, T. Nagashima, P. Newman, S. Pawson, G. Pitari, E. Rozanov, C. Schnadt and T. G. Shepherd, 2003. Uncertainties and assessments of chemistry-climate models of the stratosphere, Atmos. Chem. Phys. 3, 1-27.        [ Links ]

Birner, T., A. Dörnbrack and U. Schumann, 2002. How sharp is the tropopause at midlatitudes?. Geophys. Res. Lett. 29, 1700, doi:10.1029/2002GL015142.        [ Links ]

Birner, T., A. Dörnbrack and U. Schumann, 2003. The tropopause sharpness in the extratropics, Geophysical Research Abstracts. 5, 12486.        [ Links ]

Brühl, C., 1987. Ein effizientes Modell für globale Klima- und Luftzusammensetzungsänderungen durch menschliche Aktivitäten, Doctoral Thesis, Mainz, 1987.        [ Links ]

Brühl, C., P. J.Crutzen and J. U. Grooâ, 1998. High-latitude, summertime NO_{x} activation and seasonal ozone decline in the lower stratosphere: Model calculations based on observations by HALOE on UARS. J. Geophys. Res. 103, 3587-3597.        [ Links ]

Brühl, C. and P. J.Crutzen, 1988. Scenarios of possible changes in atmospheric temperatures and ozone concentrations due to man's activities, estimated with a one-dimensional coupled photochemical climate model, Climate Dyn. 2, 173-203.        [ Links ]

Fiebig-Wittmaack, M. and W. Börsch-Supan, 1994. Positivity preserving in difference schemes for the 2-D diffusive transport of atmospheric gases, J. Comput. Phys. 115, 524-529.        [ Links ]

Gallardo, L., A. M. Córdova, M. Rojas, J. Quintana, R. Alcafuz and I. Ramos, 2004. Stratosphere-troposphere exchange (STE) over Easter Island and Cerro Tololo stations in South America: analysis and simulation of an intensive sounding campaign. 8th Internacional Global Atmospheric Chemistry Conference, Christchurch, New Zealand, September 2004.        [ Links ]

Gidel, L. T., P. J. Crutzen and J. Fishman, 1983. A two-dimensional photochemical model of the atmosphere, J. Geophys. Res. 88, 6622-6640.        [ Links ]

Grooâ, J. U., C. Brühl and T. Peter, 1998. Impact of aircraft emissions on tropospheric and stratospheric ozone. Part I: Chemistry and 2-D model results, Atmos. Environ. 32, 3173-3184 .        [ Links ]

Ko, M. and A. Douglass, 1993. Update of model simulations for the effects of stratospheric aircraft. In: The atmospheric effects of stratospheric aircraft: A Third Program Report, 209-243. NASA Ref. Publ. 1313, Washington, D. C.        [ Links ]

Marchuk, G. I., 1975. Methods of numerical mathematics, Springer-Verlag, Berlin/Heidelberg/ New York.        [ Links ]

Marchuk G. I. and Skiba Y., 1976. Numerical calculation of the conjugate problem for a model of the thermal interaction of the atmosphere with the oceans and continents, Izvestiya, Atmospheric and Oceanic Physics. 12, 459-469.        [ Links ]

Marchuk G. I. and Skiba Y., 1992. Role of adjoint equations in estimating monthly mean air surface temperature anomalies, Atmósfera 5, 119-133.        [ Links ]

McRae, G. J., W. R. Goodin and J. H. Seinfeld, 1982. Numerical solution of the atmospheric diffusion equation for chemically reacting flows, J. Comput. Phys. 45, 1-42.        [ Links ]

Rondanelli, R., L. Gallardo, R. Garreaud, 2002. Tropospheric ozone changes at Cerro Tololo (30°S, 70°W, 2200 m) in connection with cut-off lows. J. Geophys. Res. 10.1029/2001JD001334, 03 December 2002.        [ Links ]

Schneider, T., 2004. The tropopause and the thermal stratification in the extratropics of a dry atmosphere, J. Atmos. Sci. 61, 1317-1340.        [ Links ]

Wirth, V., 2004. A dynamical mechanism for tropopause sharpening. Met. Zeitschrift, submitted.        [ Links ]