Introduction
Many kinds of airborne impurities may be observed in the atmosphere. Examples are desert or volcanic dust, water isotope species, and radioactive and chemical pollutants. The study of the atmospheric cycles of such tracers is interesting from many viewpoints: the potential danger to ecosystems represented by some chemical tracers, such as acids, is in itself a strong motivation; other tracers, such as HDO and H2 18O, are important as palaeoclimatic tracers; radioactive tracers (tritium, radon, etc) are useful for determining air mass trajectories; finally, aerosols may have an important impact on climate by modifying the radiative budget of the atmosphere. For these various reasons, a lot of effort has been devoted in recent years towards analyzing and monitoring the distributions of such tracers, although the data sets often remain sporadic. It may be timely to take advantage of the reliability now achieved by general circulation modelling of the atmosphere and try to supplement these data sets by numerical simulations of global airborne impurity cycles. This approach is particularly useful for studying such cycles in palaeoclimatic conditions where data sets are particularly scarce. For instance, the relations between the different water isotope species, HDO, H2 18O and H2 l6O, and climatic parameters such as surface air temperature and precipitation, are relatively well known under present conditions (Reference DansgaardDansgaard 1964, Reference Lorius and MerlivatLorius and Merlivat 1977), However, their variation with climate is unknown and this prevents the full use of isotope data for palaeoclimatic reconstructions. Another example concerns desert dust particles whose palaeo-deposition rates are known mainly from observations at a few stations on ice caps; the very large deposition rates observed at the last glacial maximum (18 ka BP) (Reference Petit, Briat and RoyerPetit and others 1981) are therefore difficult to interpret in the absence of a global numerical model simulating variations in remote continental source areas and variations in transport processes.
A Brief Survey of General Circulation Modelling
Atmospheric general circulation models (AGCMs) simulate the time evolution of various atmospheric fields (wind speed, temperature, surface pressure, specific humidity) discretized over the globe, through the integration of the basic equations: hydrostatic equation of motion, thermodynamic equation, mass continuity equation and water vapour transport equation. To reproduce the observed regime of atmospheric circulation, these dynamical equations must be supplemented with a modelling of the various sources and sinks which drive it. Therefore AGCMs also include parameterizations of radiative budgets, fluxes of momentum, latent and sensible heat at the surface, release of latent heat through condensation in oversaturated areas and various internal processes which operate at scales smaller than the relatively large mesh size of the model. These latter processes include three-dimensional turbulence in the boundary layer which transfers fluxes of momentum, latent and sensible heat upwards to the free atmosphere, lateral diffusion through subgrid-scale horizontal motion and cumulus convection which drives convective precipitation and redistributes momentum, heat and water vapour in the whole cloud column. In order to compute surface budgets, parameterizations of land or ice surface reservoirs of heat, water content and snow cover are generally included, while surface albedo, sea surface temperature, sea ice and glacial ice extent may be prescribed at their climatologicat mean value. These prescribed boundary conditions plus incoming solar radiation and CO2 content drive the atmospheric circulation towards mean climatic conditions and must be modified when considering seasonal variations or palaeoclimatic conditions. Integration is performed starting from realistic initial states; the first simulated days then correspond to a forecast experiment. After about ten days, deterministic predictability is lost by turbulence and long-term simulations are analyzed in terms of statistical averages, i.e. climatic averages, AGCM performances for reproducing seasonal climatic variations are quite good and may be used to perform climatic experiments or for simulations of airborne impurity cycles in the atmosphere.
Airborne Impurity Cycle Modelling
AGCMs are appropriate for studying large-scale transport of airborne impurities with residence time scales varying from days to months. Note however, that stratospheric circulation is not very realistic in most AGCMs; stratospheric pollutant studies therefore necessitate special models with higher stratospheric resolution and more elaborate treatment of stratospheric physical processes. One possible approach is to separate the pollutant transport from the basic AGCM simulation, using the simulated wind fields as input for the tracer model; the pollutant transport can be computed in Eulerian form (advection of tracer density, e.g. Reference Mahlman and MoximMahlman and Moxim (1978) or G L Russell, J A Lerner and E M Moore (private communication)) or in Lagrangian form (calculation of trajectories of individual particles, e.g. Kida (Reference Kida1983, Reference Dansgaard, Johnsen, Clausen and Gundestrup1984)). The other approach is to treat the tracer as a supplementary variable in the AGCM itself. The former approach is more economical because a variety of tracer experiments can be performed from a single AGCM simulation. However, the latter approach is more accurate because a number of important dynamical processes, for instance, convective diffusion, cannot be precisely reconstructed from low-resolution time series of AGCM basic variables. This is especially true of isotope cycles that have to be modelled following closely the basic water vapour cycle of the model. The latter approach, on which we shall concentrate here, has been utilized by Reference JoussaumeJoussaume (unpublished) and Joussaume and others (Reference Joussaume, Sadourny and Jouzel1984[a] Reference Joussaume, Rasool, Sadourny and Petit[b] Reference Joussaume, Sadourny and Jouzel[c]) on the LMD model (Reference Sadourny and LavalSadourny and Laval 1984) for desert dust and water isotope (HDO, H2 18O) cycles; a similar modelling effort for water isotopes is being carried out on the Goddard Institute for Space Studies (GISS) AGCM (Reference HansenHansen and others 1983) by J Jouzel and others (private communication). These models can, of course, be generalized for other types of pollutants.
The evolution in time of a tracer is governed by source and removal processes, which may be highly dependent on the nature and type of particles, and on large-scale transport and diffusion processes. All gaseous and particulate pollutants will be considered here as passive scalers which means that their eventual effect on climate, either as condensation nuclei or through their radiative properties, will not be included.
1 Source terms
From one pollutant to another, we may find quite different source mechanisms, varying, for example, from volcanic eruptions, fires, and dust storms to anthropogenic or extraterrestial origin (see Reference Caddie and RasoolCaddie 1973 for more details), and thus encounter a wide range of modelling problems. In each case, both input flux and source areas must be either modelled or prescribed. Note that the large mesh size considered in AGCMs may lead to a coarse representation of some source mechanisms occurring on a small scale.
1(a) Dust particles
For particles raised by winds, e.g. dust over desert areas, the surface turbulent flux can be parameterized using a bulk aerodynamic formulation similar to the one commonly used for momentum, latent and sensible heat fluxes. In that approach, fluxes depend on surface roughness through a drag coefficient which may depend on particle type, on the magnitude of surface wind velocity, and on the difference between saturation mixing ratio at the surface and particle mixing ratio at 10 m above.
In the LMD model, the same drag coefficient is used for momentum, latent and heat flux; this assumption has been extended to desert dust particles. We may consider that particle loading near the surface is limited by the available turbulent kinetic energy: in that case, the saturation mixing ratio at the surface becomes proportional to the square of surface wind velocity (Reference JoussaumeJoussaume unpublished); it can be shown that this approach agrees with Reference Gillette and MoralesGillette’s (1977) experimental law. At present, we have no real estimate of the proportionality constant; specifying this constant, however, amounts to choosing a particular mixing ratio unit since all processes involved are linear.
Realistic source regions for fine dust should be prescribed a priori including, for instance, dried-out alluvial plains. Reference JoussaumeJoussaume (unpublished), however, did not prescribe them but simply defined them as regions where instantaneous soil moisture gets lower than a given threshold value. This definition is too simple to be realistic because it deliberately ignores soil type and vegetation. It was chosen because the model is intended to be used for palaeoclimatic reconstructions where dust source regions are unknown. The implicit assumption is that source areas are distributed within, and proportional to, dry areas.
1(b) Water isotope species
Water isotope surface fluxes are computed similarly to water vapour flux, assuming isotopic equilibrium over oceans and no fractionation over land. Like most AGCMs, the LMD model uses a simplified parameterization of soil moisture with no consideration of deep soil-water content. For water isotope content, we even simplify further by considering a single mixed-layer reservoir with no variation of ground isotopic ratio with depth (Reference Joussaume, Sadourny and JouzelJoussaume and others 1984[a] Reference Joussaume, Rasool, Sadourny and Petit[b] Reference Joussaume, Sadourny and Jouzel[c]). These assumptions are crude and might not be adaptable to accurate simulation of seasonal cycles.
1(c) Other examples
Volcanic aerosol sources are easier to model as prescribed instantaneous local inputs. For radioactive gas, such as radon, the assumption of constant flux over continents seems fairly realistic (G Polian and J Sanak private communication).
2 Removal processes
Various removal processes occur in the atmosphere. Some gaseous pollutants are sensitive to chemical reactions or radioactive decay (e.g. see Reference Hidy and RasoolHidy 1973). For particulate pollutants, the dominant removal processes depend largely on size. For large particles (>10 μm) sedimentation prevails, for medium-size particles (around 1pm) rainout, in which particles act as condensation nuclei, is the main process, for Aitken nuclei (<10−2μm) coagulation of particles dominates (see Reference Hidy and RasoolHidy 1973 for more details). The spectrum of desert dust in the atmosphere lies mainly between 10−2 and 100 μm with a maximum around 10−1 μm (Reference D’Almeida and SchützD’Almeida and Schütz 1983) but only particles in the range from 10”2 to I μm will be transported far from their sources. Therefore bigger ones have been neglected and only one size range of particles (from 10−1 to 1 μm) has been modelled. Rainout processes have been parameterized with the simple assumption that dust particles and water vapour are removed in the same proportion by condensation (Reference JoussaumeJoussaume unpublished).
Removal processes for water vapour are governed by condensation processes. The parameterizations used depend on the type of cloud considered: for instance, stratiform clouds or convective clouds. The treatment of water isotope species closely follows these parameterizations with the added complication of systematic isotopic fractionation at phase changes (Reference Joussaume, Sadourny and JouzelJoussaume and others 1984[a]).
3 Diffusion processes
The two main diffusion processes are turbulent diffusion, especially vertical diffusion within the planetary boundary layer, and convective diffusion within convective clouds. In Reference JoussaumeJoussaume (unpublished), turbulent diffusion of pollutants within the planetary boundary layer is modelled exactly like water vapour diffusion, using a turbulent diffusion coefficient which is a function of static stability and vertical wind shear (Reference Sadourny and LavalSadourny and Laval 1984). The modelling of convective diffusion for pollutants is highly dependent on the convective cloud parameterization used in the basic model. The more elaborate parameterizations actually compute vertical mass fluxes within clouds (Reference Arakawa and SchubertArakawa and Schubert 1974). It is then straightforward to model the convective diffusion of particles and gases. In simpler parameterizations (Reference KuoKuo 1965, Reference Manabe and StricklerManabe and Strickler 1965), the dynamics and thermodynamics of individual air parcels remain unknown so that convective diffusion of pollutants has to be designed in a somewhat arbitrary way. The LMD model corresponds to this latter case (Reference Sadourny and LavalSadourny and Laval 1984). In Reference JoussaumeJoussaume (unpublished) total mixing of dust particles, water vapour and isotope species is assumed in a first stage within the whole cloud: the resulting high convective diffusion is compensated, at least partially, in a second stage by intensive condensation in upper cloud levels and intensive re-evaporation of precipitating water in lower cloud levels. The total mixing assumption would indeed correspond to maximum possible diffusion only for particles not removed by rainout like gaseous particles: the method should eventually be adapted for that case.
4 Transport processes
Transport of positive quantities is a critical point to model: ordinary numerical schemes used to discretize advection terms with finite difference or spectral methods hardly ensure positive values. This bias can be compensated by diffusion. The first order upstream scheme for instance, which is strongly diffusive along the flow streamlines, indeed ensures strict positive values. Higher order schemes in general do not ensure positive values at all; an ad hoc diffusion has to be added just to avoid negative values. Elaborate schemes such as corrected flux transport (CFT) techniques or total variation decreasing (TVD) schemes (see also the slope scheme of Reference Russell and LernerRussell and Lerner (1981)) have been developed, which ensure positive values with the least amount of diffusion; most of them are, however, expensive.
The transport of water isotope species is a particularly critical point: all three isotope species (H2 16O, H2 18O and HDO) must be transported numerically in a consistent way. Improperly chosen schemes, which do not ensure a minimum consistency, quickly lead to instabilities (isotopic ratios up to 108‰ have been obtained). However, simple schemes can be found which lead to relatively satisfactory results (Reference JoussaumeJoussaume unpublished).
An Example of Numerical Simulation
Preliminary simulations of desert dust and HDO, H2 18O cycles have been performed for a January climate with the LMD AGCM using a low-resolution version (32 points in longitude, 24 points in latitude, 11 layers, which corresponds to a horizontal mesh size of the order of 800 × 800 km2 around 50° and 1250 χ 450 km2 at the equator). All results presented here concern averages over the last 20 d of these runs which lasted 60 and 80 d, respectively.
1 Desert dust cycles
In order to obtain more satisfactory information about the transport of desert dust, we differentiate between six types of dust particles according to their geographical origin (North America, South America, Sahara, southern Africa, Eurasia, Australia); this separation is possible since all mechanisms involved are linear. The extent of simulated mobilization areas is realistic (Reference Joussaume, Rasool, Sadourny and PetitJoussaume and others 1984[b]) except for overestimation in south-east Asia where vegetation cover has not been included. Saharan dust plumes (Reference Joussaume, Rasool, Sadourny and PetitJoussaume and others 1984[b]) are reasonable, especially the Atlantic plume towards South America, in agreement with Reference Prospero, Glaccum and NeesProspero and others (1981). Figure 1 shows simulated Australian dust plumes at 900, 500 and 200 mbar: transport is dominated by westerly winds in the roaring forties, its northward extent being limited by the Indonesian monsoon, Australian dust removed by rainout (Fig.2) is mainly precipitated over Australia, where dust concentrations are high, and in regions of high precipitation rate such as Indonesia and the roaring forties; however, it must be pointed out that precipitations in the roaring forties are underestimated in this low-resolution version of the model.
Fig. 1 20-d average simulated Australian dust plumes at (a) 900, (b) 500 and (c) 200 mbar Shadine interval: 0.11, 0.24 × 10−1, 0.39 × 10−3, respectively.
Fig. 2 20-d average simulated Australian dust precipitations. Shading interval: 0.38 × 102.
As already noticed, dust mixing ratios and therefore dust deposition amounts are defined in our model up to a multiplicative constant which forbids any quantitative estimate of absolute deposition rates. The use of different dust types, however, allows us to estimate the relative contribution from different continents to dust deposits over ice sheets. This type of information is of obvious interest although it has to be taken cautiously since it may be particularly sensitive to the amount of diffusion included in the numerical advection schemes which modulates dust plume extents. In western Antarctica, Australian and South American dust dominate by an order of magnitude. In central Antarctica, dust originating either from Australia or South Africa is four times smaller than South American dust. Australian dust dominates by a factor of three in eastern Antarctica where South American and South African dust are also present. Greenland receives approximately equal contributions from Central America, Sahara and Eurasia; however, in this period of 20 d, the central American source area is particularly weak.
2 Water isotope cycles
A detailed description of isotope results can be found in Reference Joussaume, Sadourny and JouzelJoussaume and others (1984[a]). The geographical chart of simulated mean oxygen-18 content of precipitation (Fig.3(a)) is in good agreement with mean January observations (Fig.3(b)). Both charts show low ratios at high latitudes with the expected positive correlation with mean surface air temperature, and higher ratios in low latitudes with a negative correlation with the precipitation field. However, the model values tend to be somewhat overestimated over ice sheets: over Greenland, the simulation yields isotopic ratios around −22‰ compared with observed values down to −40‰ (Reference Dansgaard, Johnsen, Clausen and GundestrupDansgaard and others 1973); over Antarctica, the minimum simulated value is −32‰ only, instead of −40 to −50‰ in January at South Pole (Reference Aldaz and DeutschAldaz and Deutsch 1967, Reference Jouzel, Merlivat, Petit and LoriusJouzel and others 1983) and the annual averages are down to −55‰ (Reference MorganMorgan 1982). This tendency might, however, be partly due to overestimation of simulated temperatures there. We might therefore expect an improvement merely by using a higher-resolution version of the LMD model.
Fig. 3 Mean oxygen-18 content of precipitation charts, (a) 20-d average simulated field for January using the model, (b) mean January observations (J R Lawrence in preparation). Values in versus the standard mean ocean water (SMOW) reference. Isolines every 2‰ (model), and every 4‰ below −8‰ (observations). On both figures values above 2‰ are dotted and values between —8 and −4‰ are shaded.
The simulated and observed s18O-temperature relationships (Figs. 4(a) and (b)) are also in quite good agreement; the higher slope and higher correlation obtained in colder regions from our simulation agree particularly well with observations. Observations (Fig. 4(b)) also show that isotopic ratios over ice sheets tend to be lower than those given in the statistics of the International Atomic Energy Agency (IAEA); this is illustrated by Dansgaard’s line (Reference DansgaardDansgaard 1964) based on annual averages; the effect seems even stronger on January data (Reference Aldaz and DeutschAldaz and Deutsch 1967, Reference Jouzel, Merlivat, Petit and LoriusJouzel and others 1983). Our simulation shows a similar tendency although simulated temperatures are warmer.
Fig. 4 Mean oxygen-18 content of precipitation as a function of mean surface air temperature, (a) simulated fields from the model, (b) mean January observations from the IAEA network. Continuous lines correspond to the regression lines:
for atmospheric temperature TA > 15°C.
where r is the correlation coefficient. The dashed line corresponds to Reference DansgaardDansgaard’s (1964) observations, δ18O = 0.69TA−13.6. Simulated Antarctic and Greenland grid point values are identified by large dots.
The δD-δ18O relationship (Fig.5(a)) is also realistically simulated, in agreement with IAEA observations (Fig.5(b)) and with Craig’s meteoric water line (Reference CraigCraig 1961). In particular, the effects due to kinetic fractionation during re-evaporation under cloud base (enhanced dispersions and lower slopes in regions of isotopic enrichment) are well captured by the model.
Fig. 5 Mean oxygen-18 content of precipitation as a function of mean deuterium content of precipitation, (a) simulated fields (model), (b) mean January observations from the IAEA network. Continuous lines correspond to regression lines:
where r is the correlation coefficient. The dashed line corresponds to Reference CraigCraig’s (1961) observations:
Conclusion
Although the simulations done so far with desert and water isotopes are examples only, they seem to demonstrate the capability of AGCMs to reproduce global airborne impurity cycles. Such simulations, supplementing observations, should give us the opportunity to improve our knowledge of gaseous and particulate pollutant cycles with a possible feedback on our understanding of air mass trajectories. There remains, of course, a number of critical modelling problems which need further improvements: numerical advection schemes and convective parameterizations, for instance, should be modelled in a more accurate way. Quantitative improvements of models, however, require more global and more quantitative observations, especially concerning desert dust mixing ratios, deposit distributions and mobilization areas.
Acknowledgements
This work has been realised with R Sadourny (Laboratoire de Météorologie Dynamique) and J Jouzel (Laboratoire de Géochimie Isotopique) and 1 thank them particularly for their essential contributions. I am also grateful to Ρ Le Van (LMD) for his important help with the programming work, and to G Rabreau (LMD) for developing the graphics. I thank J R Lawrence (Lamont Doherty Geological Observatory) for providing us with the geographical chart of mean observed January oxygen-IE content of precipitation. This work was supported by the Programme National de la Dynamique du Climat.
