Introduction
The 18O/16O ratio is given in terms or the unitless variable δ(18O) (in ‰, and referred to hereafter as δ) and has been used extensively to study the movement of atmospheric water vapour. Oxygen isotope records from ice cores have been used previously to provide local temperature histories at polar sites (Reference Dansgaard, Johnsen, Clausen and GundestrupDansgaard and others 1973, Reference Barkov, Korotkevkh, Gordiyenko and KotlyakovBarkov and others 1977, Reference Budd and MorganBudd and Morgan 1977, Reference PatersonPaterson and others 1977, Reference Lorius, Merlivat, Jouzel and PourchetLorius and others 1979). It is known that up to seven non-temperature effects can alter the δ content of precipitation at a given site (Reference Dansgaard, Johnsen, Clausen and GundestrupDansgaard and others 1973, Reference FisherFisher 1979), but the assumption has held that the local site temperature remains as a strong(est) signal in the δ time series. This assumption rests on a view that most water vapour arriving at polar sites is transported from the warm tropics (10 to 25°N or S), where evaporation rates and sea surface temperatures are high (Reference Dansgaard, Johnsen, Clausen and GundestrupDansgaard and others 1973, Reference Lorius, Merlivat, Jouzel and PourchetMerlivat and Jouzel 1979). In the simplest form of the theory, isotopic fractionation en route follows the Rayleigh condensation process from a limited mass of vapour, and the δ(18O) of the condensate at a site with temperature T is
where Fv = y/y0 with y0 and y being the saturation mixing ratios of water vapour at the initial condensation temperature T0 and site temperature T. α, α0, αm are the fractionation coefficients (Reference MajoubeMajoube 1971) for 18O at T, T0 and (T+T0)/2 respectively. is δ for the initial vapour mass at T0 and is usually taken as (1/α0 − 1) which is close to −0.01 (i.e. −10δ. Other factors influence (Reference Lorius, Merlivat, Jouzel and PourchetMerlivat and Jouzel 1979), but for now we assume .
All the variables in Equation (1) are functions of either T or To, thus δ itself is only dependent on the site temperature T and the initial condensation temperature To (To = +20 °C for a sea surface temperature of +25.4°C, with 80% humidity and adiabatic cooling of the moist air as it rises). The argument continues that To in the tropics is very stable even under ice-age conditions, so that the β value of precipitation is principally related to the site temperature T.
Problems with the Theory
Where difficulties with the single-source theory have been encountered, the concept of local water vapour mixing with the tropical has been introduced. For example, Reference Koerner and RussellKoerner and Russell (1979) and Reference SchriberSchriber (unpublished: 54–77) need to add about 20% local Baffin Bay water to the annual precipitation in the Coburg Island area (76°N, 79°30ˈW, just east of Devon Island) to account for a δ too positive for the air temperature. Some difficulties are harder to reconcile with the simple theory. There is evidence (J England personal communication) that the regional air temperatures of north-east Ellesmere Island did not change sharply at the Wisconsin-Holocene δ transition at 10.5 ka BP. The 10.5 ‰ step in δ in cores from the Agassiz Ice Cap (Reference Fisher, Koerner, Paterson, Dansgaard, Gundestrup and ReehFisher and others 1983) in this area could not have been due to a major sudden warming at the site. The air temperature of northern Ellesmere Island seems to have remained cold until at least 8 ka BP. Uplift curves, dated moraines and driftwood (J England personal communication) suggest that ice retreat on northern Ellesmere Island began slowly in about 8 ka BP and speeded up around 6.2 ka BP.
Reference ErikssonEriksson (1965) pointed out that the simple theory is physically unrealistic and that a global model of δ should depend directly on many meteorological variables along the whole water cycle (evaporation, precipitation advective wind velocity, eddy diffusivity, temperature, relative humidity). He developed a set of equations containing all these variables, but data quality prevented him from solving them. Reference Robin and deRobin (1977) has examined and rejected the assumption that the region of major vapour input for polar sites is the warm tropics. He points out that little or no moisture injected equator-ward of 25° ever reached the polar regions. A number of authors have presented strong empirical evidence showing the importance of the distance from the closest moisture source in determining the δ value. For polar regions they have shown that the distance from open water is more influential than temperature in determining the δ value of the site (Reference Kato, Watanabe and SatowKato and others 1978, Reference KoernerKoerner 1979 Reference Punning, Vaykmyae, Kotlyakov and GordiyenkoPunning and others 1980, Reference Bromwich and WeaverBromwich and Weaver 1983).
The Model
The following model uses zonal averages of evaporation, precipitation, temperature, net meridional vapour flux, and the position of the sea-ice front as input. All sources allowed by the net vapour-flux field contribute towards the precipitation at a given site. Since the model does not include orographic effects and is weighted towards ocean conditions, it generates results that apply to sea-level, island and coastal stations. The impetus of the work is to produce a model that includes the major variables in the water vapour cycle and generates reasonable δ values that are not just dependent on the site temperature.
Consider a target latitude zone width Δxt (km) and distance xt (km) from the equator. The moisture falling at xt arrives from all the “upwind” zones. Assume for now that the source zones are closer to the equator than the target site, and let the amount of vapour over the target site xt that originates from a source zone centered at distance x (from the equator) be y(x,xt). Precipitation over the site removes some fraction k of all available vapour contributions. Then the total precipitation at xt is
The δ of precipitation due to vapour from a source at x is δ(x,xt) and the average δ(xt) of all the precipitation at the target site xt is
Now comes a key assumption. The moisture initially injected into the atmosphere at x, y0(x) is taken to be proportional to the area of the source zone A(x), the fraction of ocean in the zone C(x), and the oceanic evaporation rate E(x). Thus y0(x) = K. E.A.C., where K is a constant. The amount of vapour left at xt from a given source at x is
where f = y/y0, the depletion fraction.
For brevity and is called the source zone weight. Constants k and Κ have disappeared in the important weights in Equation (3)
The weights for the northern hemisphere are obtained from annual zonal averages and appear in Table I for present-day conditions.
The ocean-weighted evaporationrate and the ocean fraction are used to calculate the weights instead of where is the total zonal evaporation because the model is attempting to generate sea-level, island and coastal δ values. The assumption behind this is that the most precipitation at such sites originates from the oceans. This is reasonable because most of the northward net meridional vapour flux in the northern hemisphere is over oceans (Reference Peixóto and StarrPeixóto and Starr 1958: fig. 112). In practice the results are not very sensitive to whichever choice of E is used.
The depletion factor f is derived using the vapour survival time t(x), i.e. the average time that vapour at a given latitude stays aloft. During time t(x), vapour at x can travel some distance, whose meridional component is denoted λ(x), and can be thought of as the average meridional survival distance for water vapour at x. There are various estimates of t(x). We use that given by Reference JungeJunge (1963) based on t(x) = W(x)/P(x), where W is the amount of precipitable water and Ρ is the precipitation rate. We define λ(x) from the net annual zonally-averaged meridional vapour velocity V ϕ:
V ϕ has been defined by
where S Q ϕ is the net annual zonally-averaged meridional vapour flux across a latitude line at distance x. W is the precipitable water vapour and S the peripheral length of the latitude circle at x. Included in S Q ϕ are both advective and turbulent, or diffusive, vapour transfer. Since S Q ϕ includes transport at all elevations the model will not produce elevation effects in the results. The model should be compared to sea-level observations only. Estimates of S Q ϕ and W come from Reference Peixóto and StarrPeixóto and Starr (1958) and are listed in Table II along with t(x) and λ(x).
Figure 1 shows the mean annual meridional survival distance λ(x) based on data for one year (Reference Peixóto and StarrPeixóto and Starr 1958). Positive values of λ indicate northward vapour transport and negative values southward. This figure has two node points So and No, where λ = 0. More recently published data by Reference Peixóto, Oort, Street-Perrott, Beran and RatcliffePeixóto and Oort (1983), based on data for ten years, show that λ remains large and positive as far as 80 °N. Their data also show the maximum of λ in the northern hemisphere to be shifted south, and the So point to remain as it is. Under the formalization of this model precipitation at So all originates locally. In reality, the nodal point moves and precipitation at So comes from sources north and south of So.
The following approximation is used to estimate λ(x) north of So in Figure 1:
where x* is the symmetry point (= 5820 km), λ* is the maximum λ (= 1450 km), and B is half the width between the intercepts of the parabola (= 3600 km), The λ* value from Figure 1 is seen to be 900 km, but the larger value above is used (the figure of 900 km is a total zonal average). Over the oceans the northward vapour flux is larger and thus the λ* appropriate to the ocean-weighted model is larger, besides by experiment the constants listed for Equation (7) give the best results. The constants of Equation (7) result in positive values of λ between 21 and 85°N. If one could follow an amount of water vapour yo(x) originating from a zonal band centered at x > So, it would move northward and be depleted according to
If λ were a constant one would get the simple e-folding depletion relation y/yo = e−x/λ. With λ(x) given by Equation (7) one gets the depletion fraction at target xt of vapour originating at x
The relative depletions f = y/yo go into Equations (5), (3), and (1) so that δ at site xt can be calculated. The values of α of Equation (1) are changed here to apparent or effective fractionation coefficients αe (Reference ErikssonEriksson 1965, Reference Lorius, Merlivat, Jouzel and PourchetMerlivat and Jouzel 1979) by
Here α is “corrected” to take account of the unprecipitated liquid drop content ye in the air mass. ye/y is typically 0.2 (Reference ErikssonEriksson 1965) and is a constant here. Reference JungeJunge (1977) points out, however, that ye/y can vary from 0.1 in the tropics to over 0.5 in polar areas. The relative zonal-average precipitation rates can be calculated from Equations (9), (2) and (4).
Application of the Model
Figure 2 displays the f and δ(x,xt) functions for a target zone centered at 75 °N. The f curve shows the relative contributions of the source zones to precipitation at 75 °N, and demonstrates that all the source zones between 30 and 75 °N supply significant amounts of vapour to 75 °N. In particular the “local water” zone between 65 and 75 °N contributes about 25% of the total precipitation. This fraction is in agreement with the conclusions of Reference Koerner and RussellKoerner and Russell (1979) and Reference SchriberSchriber (unpublished) about the origins of precipitation at sea level for northern Baffin Bay.
Figure 3(a) shows that values of δ estimated from the model (solid line) run through the measured points (δ, temperature) from island and coastal stations (solid circle: Reference YurtseverYurtsever 1975; open circle: Reference Dansgaard, Johnsen, Clausen and GundestrupDansgaard and others 1973, coastal sites), and some coastal sites in the Canadian Arctic Islands. Comparison to marine values of δ is appropriate because of the model parameters used and the dominant meridional vapour flux being mainly marine. The points from Yurtsever’s plot that are described as “suffering” from the amount effect (Reference DansgaardDansgaard 1964) are not included in Figure 3(a). The model temperatures are those zonal averages taken from Reference SellersSellers (1967) and listed in Table I. Model values of δ vs. latitude appear in Figure 3(b). No effort is made here to calculate values of δ outside the latitude range 20 to 80°N, although this can be done by inserting different constants in Equation (7) and remembering that the vapour flux is from north to south. One can see immediately that if the λ(x) function was constant in time, local δ maxima should occur at So because λ = 0 and all the precipitation would be local. The calculated maximum at So is −1.3‰. In reality So is an annual average position, and seasonal, secular, and random movement of So would produce a blurred or scattered maximum at So. The Greenland ice-cap line δ(18O) = 0.62 T – 15.2‰ (Reference Dansgaard, Johnsen, Clausen and GundestrupDansgaard and others 1973) is plotted as a dashed line in Figure 3(a). These stations would not be expected to fall on the model line or the coastal marine line because the dδ/dT relation for them is probably determined by the orographic-adiabatic mechanism described by Reference DansgaardDansgaard and others (1964, Reference Dansgaard, Johnsen, Clausen and Gundestrup1973). Thus, the values of δ on the ice cap are determined by the δ of the vapour delivered to the coast of Greenland by the model, and then by the mechanism for adiabatic orographic depletion. Once the top of an ice cap is reached the non-adiabatic depletion mechanisms can become important again (Reference Koerner and RussellKoerner and Russell 1979).
Figure 3(c) shows that the predicted precipitation (solid line) matches the measured values (broken line. Reference SellersSellers (1967)). The curves are all normalized because the model produces relative values of P(xt). In the Queen Elizabeth Islands, Canada (latitude >72 °N), the predicted precipitation matches more precisely the regional precipitations measured by Reference KoernerKoerner (1979) (see crosses on Fig.3(c)). The predicted precipitation should be close to the actual values because the t(x) function from Reference JungeJunge (1963) is calculated using zonal averages of precipitation. The agreement is thus a check of internal consistency.
Seasonal δ Variations
The seasonal latitude cycle in temperature T(x), evaporation E(x) and e-folding length λ(x) will produce a seasonal δ cycle at any given target site. For a first look at the predicted seasonal δ cycle, the annual average λ(x) function is used along with estimates of seasonal zonal averages of T, which are the seasonal temperature amplitudes in Reference LambLamb (1972: 145) weighted towards the oceans. The annual average E field is used except that E = 0 over sea ice. Since most of the measurements of the seasonal δ cycle north of 50 °N are from stations on either side of Baffin Bay, we decided to use the model in this context by using the seasonal sea-ice cycle from Baffin Bay. In January the sea ice reaches the southern tip of Greenland (60°N), and in July the bay is free of ice with the nominal northern limit of open Baffin water at 75 °N. Thus, there is a 15° latitude sea-ice cycle. Also, we assume that sites north of 75 °N in the Queen Elizabeth Islands and northern Greenland receive no moisture from sources between 75 to about 83 °N, this region being mostly land.
Figure 4 summarizes the measured seasonal amplitude A = (δsummer – δwinter)/2 for Greenland and Canadian Arctic ice cap stations, coastal Greenland stations and some Antarctic sites. The dashed lines are linear regressions through the northern ice-cap and coastal stations. The solid line is the model prediction. The model line follows the coastal station line (lower dashed) until about 65°N, and then rises to meet the ice-cap line near 75°N where A reaches a maximum for both model and ice cap. The contribution of the sea-ice cycles to A can be seen immediately by comparing coastal stations at either limit of the sea-ice cycle. At Gronnedal (61 °N), which feels little or no sea-ice effect, A = 3.3‰ and at 75 °N, which gets the full effect, A = 6.4‰. Thus the maximum contribution to A of a 15° latitude march in sea ice is about 3.1‰, about half of the total A at 75 °N. With the rough seasonal parameters used here the model predicts a 4.8‰ contribution to A from the sea-ice cycle and 4‰ from the other seasonal variables. The rise of A in the model above the measured coastal line north of 65 °N is possibly because the seasonal values of E used are not realistic north of 60°N; also one cannot hope to reproduce much more than general tendencies from a zonal model even if it is “tuned” to a particular area. The “odd” small measured values of A for the ice-cap sites north of 80° could be related to the presence of the No mode at 77 °N. Sites north of No receive much of their moisture from the Arctic Ocean.
It is interesting to note that for the Antarctic Reference Robin, de, Johnsen, Robin and deRobin and Johnsen (1983) have drawn attention to the fact that A at the coastal site of Syowa (latitude 69 °S) is close to that of Mizuho camp (latitude 70°45ˈS) (Reference Kato, Watanabe and SatowKato and others 1978) and Amundsen-Scott station (latitude 90 °S) (Reference Aldaz and DeutschAldaz and Deutsch 1967). The coast in this case is the southern limit of the sea-ice cycle. Robin and Johnsen go on to suggest that the biggest contribution to A is caused by this cycle. We find this assertion too strong for the northern hemisphere. At 75°N only 50% of coastal A is due to the sea-ice cycle. The values of A in the Antarctic are plotted in Figure 5 with the Syowa point placed at 75 °S, which is the poleward limit of the northern sea-ice cycle. The A-latitude relationship poleward of 80° might be different for the two hemispheres because in the north there are moisture sources north of 80°, but none in the Antarctic south of 70°.
Model Applied to Hypothetical Ice-Age Conditions
In order to try out the model on some hypothetical conditions 18 ka ago, we assume that the annual average λ(x) is as today. T(x) is derived by altering today’s temperature averages according to the CLIMAP 18 ka BP temperature for the Atlantic Ocean (CLIMAP Project Members 1976). In the high Arctic the 18 ka BP temperature is obtained from the climatic δ shift of the Devon Island ice core (Reference PatersonPaterson and others 1977). This 7‰ shift suggests a temperature 11 °C cooler than today. Although this use of δ to estimate the 18 ka BP temperature is circular, we note that the 11 °C fits in as a reasonable northward extrapolation of the CLIMAP ocean temperature changes. Also, the values of δ in the model at high Arctic target sites are not markedly sensitive to the temperature at the target site. These 18 ka BP temperatures are listed in Table I, In addition, we assume an 18 ka BP sea-ice cycle having the January ice front at 55°N and July at 65°N. The annual average position is about 10° farther South than today. The calculated mean annual values of δ are shown in Figure 3(b) plotted against latitude, and the inset shows the (δϕk-δ18k) difference. The δ differences are about the right amount. Table III lists this difference for the ice cores drilled north of 60 °N. The δ difference for the Barnes Ice Cap (Reference Hooke and ClausenHooke and Clausen 1982) is complicated by extreme changes in the elevations of the ice-flow line, present-day melt on the surface and the non-continuous δ profile available. Recent flow-model calculations suggest that during the late Wisconsin the flow tines through the Barnes area originated 2400 m a.s.1. (N Reeh personal communication). This underlines the caution involved in comparing the predictions of a zonally-averaged sea-level model to the differences measured at ice-cap sites at high elevations where so many other effects related to the ice cap can influence the values of δ. This rough model only indicates the rough order of magnitude of the δ difference and suggests that (δø–δ18k) has some inherent latitude variance.
Deuterium-Oxygen 18 Relationship
Up to this point we used Equation (1) which gives an approximation to the precipitation S produced by condensation from a limited quantity of vapour from a given source zone. Also we have assumed that the initial vapour over the source zone to be a constant −10‰. In order to extend this model and examine the δ18Ο – δD relationship Equation (1) is replaced by the differential equation developed by Reference Lorius, Merlivat, Jouzel and PourchetMerlivat and Jouzel (1979).
where, as before, y is the mixing ratio of water vapour, ye is the mixing ratio of liquid droplets in the clouds, and α is the fractionation coefficient for 18O or D.
Their expression for the initial vapour δ is
where αc is the isotopic fractionation coefficient (18O or D) at the sea-surface temperature and h is the relative humidity in the air mass over the ocean, k is the variable related to kinetic (non-equilibrium) fractionation effects in the air-ocean boundary layer. For 18O
and
(Reference Merlivat and JouzelMerlivat and Jouzel 1979: fig.2). In addition, they show that for deuterium D
They define a rough regime as having wind speeds ≥7 ms−1 at 10 m above the surface. Although they quote Reference Eriksson and BolinEriksson and Bolin (1965) as stating that 95% of the world’s oceans are smooth, they calculate an ocean-weighted value for k18O of 4.8‰. The Present model was run with two k18O extremes of 7‰ and 4‰.
The independent variable in the integration is x, the distance from the equator, and Equation (11) is used in the form
T being the average zonal temperature and ye/y as before is kept at a constant 0.2. dy/ydx comes from Equations (7) and (8). From the region north of So integration starts at SO and proceeds northward using empirical data from Table I. Figure 5(a) shows the (δ18O, δD) points from Reference DansgaardDansgaard (1964: fig.10) (dots) along with points from time averages adapted from ice-core data, Byrd core (Reference Epstein, Sharp and GowEpstein and others 1970) (squares), and Dome C (Reference Lorius, Merlivat, Jouzel and PourchetLorius and others 1979) (ellipses). Reference DansgaardDansgaard’s (1964) least squares line is the dashed line δ(D) = 8.1 δ(18O) + 11. The model points all fall very nearly on their solid least squares line (Fig.5(a)) δ(D) = 8.37 δ(18O) + 11.9. The deuterium excess is defined by d = δ(D) – 8 δ(18O).
Reference Jouzel, Merlivat and LoriusJouzel and others (1982) examined the important variables influencing d in the ice core recovered from Dome C, Antarctica. There is a decrease of 4.5‰, in d from Holocene to ice age. They concluded that (i) the d excess is conserved from source to target, (ii) d is mainly determined by the relative humidity h of the air over the source ocean, and (iii) the d shift of 4.5‰ from present to ice age is thus due to a change in relative humidity from 80 to 90%.
Using the modern atmospheric and oceanic data for the southern hemisphere and an annual average sea-ice front at 65 °S (Ackley 1981), the present multi-source zonal model calculates the deuterium excess shown in Figure 5(b) for 80 to 90% relative humidity and for smooth and rough source oceans. Since both h and k are probably latitude variables the actual d values are likely to be some average of the curves. The sparse “measured” values of d for the southern hemisphere are plotted. Since the zonal model is weighted to the oceans and relevant to sea-level, the model values of d are not really valid for latitudes south of 70°S. It seems from Figure 5(b) that d presently has a latitude dependence and, furthermore, at any given site d can be shifted as much as 8‰ by altering k and h within reasonable limits.
Speculative runs of the model for the southern hemisphere were made assuming ice-age (18 ka BP) air and sea temperatures altered by amounts suggested by CLIMAP Project Members (1976: fig.3) and Reference Robin, de, Johnsen, Robin and deRobin and Johnsen (1983: chapter 6.3) and taking an average ice-age sea front at 55 °S. The calculated d values suggest that half of the 4.5 ‰ d change at Dome C could be related to the latitude dependence of d. Because of other evidence presented by Reference Jouzel, Merlivat and LoriusJouzel and others (1982) (i.e. higher ice-age salt concentrations and atmospheric model predictions of enhanced storminess), it seems that k, h and their latitude dependence are important in explaining the d shift. Assigning all the d shift to an h change seems to be arbitrary. Of course, a zonal model can only suggest what variables are important in affecting d at any given site. For example, Figure 5(b) shows a change in d in the Byrd station core of about ‰. The zonal model cannot explain such differences from one longitude to another. Finally, the calculated (δøk−δ18k) differences on the coast of Antarctica are about 4‰, or about half those calculated for northern ice-cap sites. This is in agreement with the ice-core data.
Further Work
In principle the theory of this simple zonally-averaged model could be carried to a latitude-longitude variable model with resolution iimited by the resolution of the input data. One could then use the geographical array of ice-core 5 time series stretching back many millenia to validate and test proposed atmospheric models that produced vapour trajectories and the latitude-longitude functions of evaporation E, precipitation P, moisture flux Q, precipitable water W and temperature T. Thus, if one has lost the target site δ as a simple function of the site temperature, one has gained through the geographical array of δ time series a powerful and independent set of data to test the state of the whole atmospheric flow pattern. The job is now bigger but even more rewarding.