Numerical Transport

1. Numerical Transport Anne Douglass Code 613.3 Atmospheric Chemistry and Dynamics Branch NASA Goddard Space Flight Center

2. Numerical Transport This talk is planned to convey two things: the basic performance requirements for a useful numerical transport scheme; the history of the development of these requirements, and the importance of observations to development.

3. Start with the continuity equation: Separate the solution of the advective and photochemical parts. Solver for X* using the advection operator and replace in the continuity equation: X(t) - X(0)= [X* - X(0)] t t “Process-splitting” is used in virtually all multi-dimensional atmospheric models. McRea, G. J., W. R. Goodin, and J. H. Seinfeld, Numerical solution of the atmospheric diffusion equation for chemically reacting flows, J. Compu. Phys., 45, 1-42, 1982. v ·  X + P - L

4. Once we have figured out how to obtain a solution for the advective part, we can use any of several approaches to obtaining the chemical contribution to the tendency. **all of the coupling among species is in the photochemical part** Start looking for a good solution to the advective part with a one dimensional problem: It seems logical to notice which direction the wind is blowing.

5.  v  v t = 0  v t = t t = 10 t This is “the answer” for this particular case. The square wave simply translates with velocity v. Think about transporting a square wave in the x direction with a constant velocity v.

6. We are working with a grid point model - that forces us to abandon the concept of simple translation and put this onto a grid. start out the same way - move the material 0.1 grid box to the right. remember this is what it is supposed to look like! so 90% of the material stays in the box 1, and 10% moves into the box to the right (box 2); notice that 10% of the material in box 2 will move into box 3 on the next time step ).01 of the total but it should be zero. “upwind” simply mixes the material uniformly in its grid box. box 1 box 2 box 3

7. After only one time step, we can see that such a simple scheme is not good enough. • ISSUES • The maximum drops. • The material travels too far (diffusive) CONSEQUENCES (these are modern consequences to an old problem . . .) • poor representation of transport barriers and sharp gradients? • unrealistic age spectrum and age-of-air

8. “exact” solution After many timesteps problems are obvious. Material should have moved totally from its original location, and it has not. The leading edge should be one box away, not several. palest green (lowest value in box 1) is after the 10th time step. • this example is a disaster (on purpose). I could make it better if I • increase the number of grid points • decrease the time step • *but I can not fix all the problems without changing the algorithm!*

9. Scheme adjusting means developing “flux limiters” to account for the shape of the tracer distribution within the grid box Smolarkiewicz, P. K., A fully multidimensional positive definite advection algorithm with small implicit diffusion, J. Comp. Phys.,54, 325-362, 1984. Zalasak, S., T., Fully multidimensional flux-corrected transport algorithms for fluids, J. Comp. Phys., 31, 335-362, 1979. van Leer, B., Towards the ultimate conservative difference scheme, II: Monotonicity and conservation combined in a second-order scheme, J. Comput. Phys., 14, 361-370. 1974. (and 4 more papers with the same first half title)

10. The Prather scheme is the advection algorithm that most exactly accounts for the distribution of tracer within a box: Prather, M. J., Numerical advection by conservation of second-order moments, J. Geophys. Res., 91, 6671-6681, 1986. Features: non-diffusive mass-conserving stable to Courant limit (vt/ x < 1) Disadvantage: Computational requirements (for a 3D model, nine 3D arrays for each constituent) “slopes” first order moments (overshoot and undershoot)! start with a step function upwind (very diffusive) second order moments with limiter (no undershoot) second order moments (small overshoot and undershoot)

11. Prather’s paper is unique for a few reasons: 1. published in JGR (not computational fluid dynamics) 2. actually contains code (very unusual for that era to hand out the code) “A comparison of the capabilities of different algorithms shows that each scheme has its advantages and that no one method is most advantageous under all conditions.”

12. And there were other schemes which paid no attention to the flow (did not take into account the direction of the wind): Mahlman and Moxim, Tracer simulation using a global general circulation model: results from a midlatitude instantaneous source experiment, J. Atmos. Sci., 35, 1978. 2nd order differences in the horizontal and 4th order differences in the vertical spectral transport (these schemes can produce negative mixing ratios near sharp gradients (tropopause H2O, “transport barriers” between the subtropics and middle latitudes or middle latitudes and polar vortices) “semi-Lagrangian” schemes solve both undershoots and undershoots but require attention to mass conservation Williamson and Rasch, Two-dimensional semi-lagrangian transport with shape-preserving interpolation, Mon. Wea. Rev., 117, 102-129, 1989.

13. A review paper* concluded that “. . . straightforward approaches such as spatial fourth-order finite differences with leapfrog time differences provided a cost-effective numerical algorithm. . .” First 3D CTM results from GSFC used spectral transport + Forester filter (J. Comput. Phys., 1977) + filling Rood et al., The use of assimilated stratospheric data in constituent transport calculations, J. Atmos. Sci.,46, 687-701, 1989. Kaye et al., Three dimensional simulation of HCl and HF during AASE, Geophys. Res. Lett., 17, 529-532, 1990. *Rood, R. Numerical transport algorithms and their role in atmospheric transport and chemistry models, Rev. Geophys., 25, 71 -100.

14. And there were some pretty good results from these simulations. The phrase “pretty goodresults” for a simulation only makes sense in the context of the observations available for comparison. “Kahlman filtered” CTM O3 30 hPa 2/25/79 LIMS O3 30 hPa 2/25/79 CTM O3 30 hPa 2/25/79 The “Kahlman filtered” model output looks more like LIMS data than the unfiltered output. The “pretty good” LIMS data cannot tell us if the model has problems, let alone where they come from (numerics or input winds).

15. 30N 30 hPa 45N 30 hPa 30N 100 hPa 45N 100 hPa 60N 30 hPa 75N 30 hPa 60N 100 hPa 75N 100 hPa Evolution of zonal mean ozone calculated using winds from an early data meteorological assimilation system solid lines: CTM zonal mean dashed lines: two closest LIMS points

16. 1990 Scientific Symposium “Global Tracer Transport Models” (Bermuda Biological Station for Research) Results from two dimensional models (the “assessment workhorses” of that era were nearly identical for four different transport schemes (upwind, Prather, van Leer, Smolarkiewicz) because the explicit diffusion (Kyy) used to parameterize the effects of planetary waves masked differences among transport schemes. Initial (good) results with 3D models (ours and others) and comparisons of scheme performance in 2D models contribute to the idea that “the transport schemes are good enough” and not the main sources of error in atmospheric models.

17. ‘Better’ data (spatial resolution, accuracy) from aircraft (AAOE, AASE I and II) and satellites (UARS) led to stringent testing of the CTM’s (both numerical transport and met fields contribute to transport errors). The Ozone Hole Discovery focused models on polar processes; transport schemes in spherical coordinates and the upwind grid schemes (van Leer, Prather) were stressed by the need to represent “cross polar flow”. ***Aircraft measurements show that long-lived species produce compact relationships an are highly correlated with each other. (Plumb, R. A., and M. K. W. Ko, Interrelationships between mixing ratios of long-lived stratospheric constituents, J. Geophys. Res., 97, 10145-10156, 1992. ) Some schemes cannot produce or maintain constituent correlations: Spectral transport + filters + filling 4th order differences + filling

18. We implemented the van Leer scheme in our 3D CTM process-split the transport (x direction, then y, then z) used a smaller time step at high latitudes (box sizes decrease at high latitude) to satisfy the Courant condition calculated the vertical velocity from mass continuity (assimilated winds . . . ) Allen et al., Application of a monotonic upstream-biased transport scheme to three-dimensional constituent transport calculations, Mon. Wea. Rev., 119, 1991.

19. develop the concept of a mixed “polar cap” to handle the cross polar flow problem

20. Solution 1 rev 5 rev 7 rev 8 rev Van Leer (10 revolutions) Van Leer Shape Preference (2X) 0.5 0.4 0.3 0.2 0.1 0 0.5 0.4 0.3 0.2 0.1 0 Solution 5X (0.82) 2X (0.63) 1 X (0.42) CONCENTRATION CONCENTRATION 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 LOCATION LOCATION Grid Resolution Matters! Left - keep that triangle more faithfully as grid resolution increases Right - Once the scheme has produced the shape it resolves, there is less diffusion for subsequent revolutions. PRACTICAL LESSON: aircraft exhaust, age of air place more stringent requirements on the transport scheme (and its implementation) than chlorine source gases.

21. Rood et al. Three dimensional Simulations of Wintertime Ozone Variability in the Lower Stratosphere, J. Geophys. Res., 96, 5055-5071, 1991. Bold: TOMS Data Line: Model Simulation

22. AND we produced and maintained realistic correlations among long-lived constituents. The End Not Quite the End This experience, and more and more observations to compare with simulated fields, led to articulation of new requirements for numerical transport schemes.

23. From Lin and Rood [1996] • conserve mass without a posterior restoration (implications for “numerical continuity” and use of assimilated winds • compute fluxes based on the subgrid distribution in the upwind direction • generate no new maxima or minima • preserve tracer correlations • be computationally efficient in spherical geometry (allow vt > x) • we would like to ascribe transport errors to errors in (assimilated) winds, not to the numerical scheme

24. Other things to keep in mind Convergence - Michael Prather has been doing some tests with his scheme vs. Lin and Rood. Do simulations of SF-6 (the age of air calculation) converge as grid resolution is increased? (yes - but there are differences between ‘practical’ grid resolution and the very high resolution calculations) Do the two transport schemes converge to the same result? (probably- but differences between the schemes at practical resolution (2 x 2.5) are larger than the difference between the Prather scheme at 2 x 2.5 and higher resolution - that is, the Prather scheme converges faster).

25. HALOE CH4 SEPT OCT NOV On-line 2 x 2.5 4 x 5 RESOLUTION!! CTM runs - SAME WINDS and methane loss at 4x5 the CTM mixes too much! HALOE shows separation (limited mixing) CTM ok at 2 x 2.5

26. How important is it that the CTM (or coupled chemistry climate model) reproduce the separation between the vortex and middle latitudes? (obviously a rhetorical question)

27. Lower stratospheric chlorine mixing ratio depends on both the horizontal resolution AND characteristics of the input winds GCM winds, 2 x 2.5 resolution GCM winds, 4 x 5 resolution Oct Cly (ppbv) (noisy) DAS winds, 4 x 5 resolution

28. Polar ozone evolution follows evolution of Cly. Predicted recovery is faster for DAS simulation and for GCM simulation with 4 x 5 transport calculation. Minimum 62 - 90 S (monthly zonal mean)

29. One last thought The world’s best transport scheme + the world’s best chemical scheme do not necessarily give good results Global Modeling Initiative Lin and Rood transport SMVGEAR Solver No family approximations Cl 0 at night - so at the terminator there is always a square wave. Largest transport errors, and flux limiters invoked separately for each constituent. RESULT (after 15 years of integration) - non-homogeneous distrubution of Cly (and Bry) in the upper stratosphere. (probably issues with total reactive nitrogen (NO and NO2) Family approximation eliminates the problem b/c it eliminates the square wave.

30. Conclusions • several schemes (including the Lin and Rood scheme) meet the performance criteria articulated by Lin and Rood (1996) • requirements for numerical transport performance and accuracy have been driven by observations • quality of data available for comparison with simulations • LIMS (daily maps of fourier coefficients, 4 degree latitude bands) • high resolution aircraft data (constituent correlations) • UARS (transport barriers) • SF-6 and other measures of “age” • quality of meteorological data used in the CTMs

31. Suggestion • All users of atmospheric models should devise some experiments to “exercise” the (blank!!!) model component!! • Develop personal insight into the capability of scheme • Develop a concept of errors and their importance to whatever application you have in mind • All users of atmospheric models should devise some experiments to “exercise” the numerical transport • Develop personal insight into the capability of scheme • Develop a concept of errors and their importance to whatever application you have in mind