Introduction
Stable-isotope values δ(18O), δ(D) and/or the deuterium excess d = δ(D) – 8δ(18O) in precipitation, are useful indicators of various parameters along the water cycle, e.g. Dansgaard, Reference Dansgaard1964, Reference Dansgaard, Johnsen, Clausen and Gundestrup1973; Fisher and Alt, Reference Fisher and Alt1985; Joussaume, 1985; Jouzel and others, Reference Jouzel, Russell, Souzzo, Koster, White and Broecker1987, Johnsen and others, Reference Johnsen, Dansgaard and White1989.
The model described here is empirical using zonal averages of E evaporation, Τ near-surface air temperature, Tc sea temperature, meridional water-vapour flux, h near-surface relative humidity, wind speed, sea-ice cover, Si supersaturation in the clouds, and Ts the temperature for shifting from vapour-water phase changes to vapour-ice ones. The model, like its progenitor (Fisher and Alt, Reference Fisher and Alt1985) is designed to calculate zonal averages of sea-level δ(18O), δ(D), and P, the precipitation rate. Global solutions are input to regional solutions along assumed vapour trajectories, up ice caps, over ice shelves, etc. The elevation, air temperature, and precipitation are inputs along such trajectories, and δ(18O) and d are outputs.
The main sources of the global empirical data used as input are: Oort (1983) for (T,Tc , vapour flux, wind speeds, h), IAEA/WMO (1969, 1970, 1971, 1973, 1975, 1979, 1981, 1983, 1986) for (T, h, δ(18O), δ(D)). There are supplementary sources for E from Sellers (1967), Peixoto and Oort (1983), Obcrhubcr (1988), for sea-ice cover from Jacka and others (1985), Claire and others (1987), for zonal precipitation from Jaeger (1983), and for Si(T) from Jouzel and Merlivat (1984). The database for this version of the global model is much more extensive and accurate than the first version in Fisher and Alt (1985).
The Model
Total water balance
Over a given target zone of width Δxt (km) some distance x t from the Equator, the water vapour aloft comes from all the “up-wind” source zones. The amount of vapour over xt that originates from a source centred at x (km from the Equator) is denoted y(x,xt). The initial moisture amount originating at x is denoted y1(x) and is assumed to be proportional to the area of the source zone A(x), the ocean fraction C(x) and the evaporation rate E(x). Thus, y1(x) = K.E.A.C, where K is a constant. The amount of vapour left at xt that started at a given source x is:
where ƒ = y/y1 is the depletion fraction.
The total amount of water over xt is the precipitable water W(xt):
where x0 is the farthest contributing source.
The δ of precipitation due to vapour from a given source at x is denoted δ(x,xt ) and the average δ(xt ) of all precipitation at the target site is:
Assume that a “survival distance” function λ(x) exists such that an initial moisture input y1 from a strip at x is depleted according to:
If λ were constant, the depletion fraction f at site xt would just be ƒ = y(x,xt)/y1(x) = exp(–(xt – x)/λ). The actual empirical λ(x) is estimated later. Integration of Equation (4) gives:
The precipitation P(xt) at all possible equispaced target sites x0 < xt < max(xt) is proportional to:
The proportionality constant can be evaluated by requiring that the integrated precipitation equal the integrated evaporation from x0 to max(xt ). Using expressions (5) and (6), it can be shown that survival distance can be given by:
where P(x*) and λ(x*) are known values at some x*.
Although Equation (7) shows there is an analytical relation between λ(x), E(x) and P(x), λ(x) is calculated for the zonal/global part of the model from measured vapour fluxes. The model then calculates W(x) and P(x) and compares them to measured values. The reason for doing this is the lack of accuracy and consistency in the measured E and Ρ fields at mid- and high-latitudes. This causes serious errors in λ because the difference terms in Equation (7) containing E and Ρ get smaller than the present errors in the global data fields for higher latitudes.
Equation (7), however, is used to calculate the λ function for the regional solutions beyond coastlines, because then E is zero, and Equation (7) is just:
where x c is the coast or where E goes to zero, P(xc) is the precipitation out at sea before the coast is reached, λ(xc) is λ at the coast from the global solution and P(x) is the measured regional precipitation along a line perpendicular to the precipitation isopleths, also assumed to be a vapour trajectory.
The empirical survival function
In the first version of this model (Fisher and Alt, Reference Fisher and Alt1985), λ(x) was calculated using:
where q is the mixing ratio (g water vapour/kg air), v is the average meridional velocity of the atmosphere (m s−1), W is the precipitable water in meters, and Ρ is the precipitation rate in m s−1. Zonal averages are denoted by brackets “[ ]”, time averages by a bar “—” and vertical averages by a hat “[Inline 1]”. A slightly different scheme is used here to estimate λ(x) but the result is much the same as from Equation (9).
This time the vapour flux is divided into three categories (Peixoto and Oort, Reference Oort1983), mean meridional, MM, transient eddy, TE, and standing eddy, SE.
where q’ and v’ are deviations from time means at a given point and q* and v* are deviations from zonal means at a given time. The total vapour flux [[Inline 2]] is their sum and the vapour velocity for each category is obtained by dividing MM, TE, and SE by [[Inline 3]], σt(q’) and σs(q*), respectively, where σt and σs are time and space variances. Then, the three velocities are weighted by the amount of flux carried by each. This final vapour velocity is multiplied by the survival time [W]/[P] as in Equation (9) to give a survival distance function.
This function is shown in Figure la for annual average conditions. When λ is positive (negative), it means the net vapour transport is north (south) and a zero at x means the vapour originating at x is all precipitated out at x. A closed water cycle goes between zeroes of λ or from a zero to a pole. For example, Figure la shows a northward moving water cycle between 20° and 90°Ν and a southward moving cycle from 20° to about 3°N. All evaporated water originating in a given cycle must precipitate out there.
Now, λ(x), as calculated above, is determined only to a “constant” because the quantities such as [[Inline 4]] are not simple products but rather covariances which contain a hidden correlation coefficient between q and v. Also is only a guide to the survival time. Here, the “raw” λ(x) function must be multiplied by a constant F found by appealing to the fact that Equations (1), (5), and (6) require that the correct λ(x) function must produce accurate simulations of the zonal averages of the precipitable water W and precipitation P. Experiments on the annual average λ(x) show that F = 2.4 gives the best fit between model W and Ρ and the measured fields (see Fig. 2). This value of F is used throughout this work. No doubt f is a function of latitude but here the observed Ρ and W fields are matched assuming that it is a constant.
Fig. 1. Zonal averages of annual input data fields versus latitude, 5° strips, (a) λ survival distance (km) for water vapour. Positive values indicate northward vapour transport, negative values southward. λ = 0 implies that all the water evaporated in that zone falls out in that zone, (b) Τ near-surface air temperature, (c) E evaporation rate, cm a−1, (d) h near-surface relative humidity, (e) k 18 kinetic evaporation fractionation factor; for 18O, k = 7‰ is for smooth oceans, 4‰ for rough oceans with wind speeds >7 m s−1 at 10 m above the surface, (f) Tc - T; difference between sea temperature Tc and air temperature T. (g) C per cent of the zone that is ocean. The arrows indicate annual sea-ice fronts.
Fig. 2. Zonal averages of annual precipitation rates, cm a−1, (a) Measured precipitation rates, (b) Modelled precipitation rates.
Isotopic balance
The differential equation used to follow δ(= δ(x,xt))(18O, or D) of precipitation for vapour from each contributing source zone is as in Fisher and Alt (1985):
where y is the mixing ratio due to a given source, and ye/y is the ratio of liquid drop content in the air mass. ye/y is kept as a constant, 0.2 after Eriksson (1965) and Junge (1977). The equilibrium-fractionation coefficients for the two species are from Majoube (1971a) for vapour to water and from Majoube (1971b), and Majoube and Nief (1967) for vapour to ice. As in Fisher and Alt (1985), dα/dx = (dα/dT)(dT/dx) and (1/y){dy/dx) comes from Equation (4).
The expression for the δ of a given source’s initial vapour comes from Merlivat and Jouzel (1979):
where T c is the source’s sea temperature, h is the relative humidity over the source ocean’s saturated boundary layer, and k is the parameter related to kinetic fractionation effects (non-equilibrium) in the air-ocean boundary. For smooth oceans k O18 = 7‰ and for rough oceans k O18 = 4‰ also k D/k O18 = 0.88. Merlivat and Jouzel defined a rough ocean as having wind speeds >7 m s−1 at 10 m above the surface.
Equation (10) is integrated from each source to each target site using the empirical data fields. Figure 1 shows the data fields used for the annual average runs and Figure 5 for the two extreme months of August and February.
Supersaturation, Si and T s
For air temperatures Τ > T s, the equilibrium-fractionation coefficients are used in Equation (10). For Τ < T s, however, α in Equation (10) must be replaced by ααk where αk is the kinetic coefficient (Jouzel and Merlivat, Reference Jouzel and Merlivat1984) and the water vapour in the clouds is assumed to sublimate directly to ice crystals instead of condensing to water droplets. For Τ < T s, the clouds are thought to be supersaturated with respect to vapour to ice phase changes. Jouzel and Merlivat (1984) developed an expression for αk and Fisher (in press) made slight improvements to it. Fisher’s version is used here:
where Si(T) = e(T)/ei(T). Τ is the ambient air temperature in the cloud, e is the partial pressure of water vapour in the cloud at T, and ei is the saturation partial pressure of vapour over ice at T (D/D’) is the ratio of molecular diffusion coefficients, [Inline 5] and [Inline 6]. α is the equilibrium-fractionation coefficient at T.
ε is close to 1 and ϒ is very close to 1. They are both weak functions of Τ and Si, In the original derivation of Equation (12), Jouzel and Merlivat (1984) did not directly include a coupled thermal equation, whereas the present version does. In practice, the results obtained for δ and d are very close to those of Jouzel and Merlivat with a tendency for the above versions of Equation (12) to give less “unstable” values for d. All of the model results presented here use T s = −10 deg from Jouzel and Merlivat (1984). The model runs reported use one of their preferred Si functions (denoted S4) Si(T) = 0.99 – 0.006T, where Τ is the condensation temperature in the clouds.
Since, for most of the global sea-level model runs, Τ > −10 the Si function is not too important but once a coast is crossed in a regional solution and Τ ≪ −10, Si becomes very important in determining d. The model δs are stable over quite a wide range of Si functions but the model ds are extremely sensitive to which Si function is used. The present model shares this property with the “single-source” model of Jouzel and Merlivat (1984) and Johnsen and others (1989). In fact, the calculated d functions for the various suggested Si functions of Jouzel and Merlivat (1984) resemble closely those authors’ versions.
The whole question of the Si functions’ effect on d is examined in detail elsewhere (Fisher, in press).
It should be pointed out that model ds are also rather sensitive to Ts. Other authors use different values, e.g. Johnsen and others (1989) use T s = −5°C. Observations of Arctic clouds by Curry and others (1989) gave −33<TS<−13°C and work reported by Mason (1975) suggested that in very clean air Ts could be about −26°C. If T s is related to the air’s microparticle loading, then considerable care should be exercised in interpretation of d time series in relation to microparticle time series in ice cores.
Global Solutions
Annual conditions
Using the empirical input functions shown in Figure 1, F = 2.4, T s = −10, and the S4 Si function, the model results are given in Figure 2b for precipitation P, Figure 3a for δ(18O), and Figure 3b for d. The dots are measured values from IAEA reports (1969, 1970, 1971, 1973, 1975, 1979, 1981, 1983, 1986) for coastal and island sites. A zonal model of course can only be expected to run through the swarm of data points. The model reproduces quite well the zonal averages of P, δ(18O) and d including the low δs around the Equator in the tropical convergence zone. In this zone between 16°S and 20°N the more negative δs go with higher precipitation, compare Figures 3a and 2b. Thus the model simulates the “amount effect”. Figure 4 plots measured and model δs against site temperatures for both hemispheres.
The model simulates the main features of the global sea-level annual d distribution with maxima at 35°N and 35°S and minima at 75°N and 70°S. The minima are produced in part by the contributions from the “local” cold oceans particularly in the more marine Southern Hemisphere.
Fig. 3. Annual average δs and ds. (a) Annual averages of δ(18O) versus latitude. The dots are measured values and the line is from the model. The model gives values near −2% where λ = 0. (b) Annual averages of the deuterium excess d = δ(D) - 8δ(18O). The dots are from measured stable isotopes, the line is given by the model using the S4 Si function with Ts = −10°.
Fig. 4. Annual average δ(18O) data (solid dots) and model δs (circles) versus site temperatures, for sea-level coastal and island stations. The hemispheres are displayed separately.
Fig. 5. Zonal averages of August (dashed) and February (solid) input data field for 5° strips, (a) λ survival distance for water vapour, (b) Τ near-surface air temperature, (c) E evaporation rate, cm year−1, (d) h near-surface relative humidity, (e) [Inline 7] kinetic evaporation fractionation factor for l8O. [Inline 8] for smooth oceans, 4‰ for rough oceans Owith wind speeds above 7 m s−1. Values between extremes of k are estimated as a linear function of wind speed, (f) Tc - Τ difference between sea temperatures and air temperatures, (g) C per cent of the zone that is ocean. Arrows indicate latitudes of sea-ice fronts.
Seasonal results from the global zonal model
The two extreme months for δ are August and February. The empirical data fields for these months (Fig. 5) were run with the same F, T s, and Si values used in the annual runs.
The model zonal precipitations (Fig. 6b) compare well with the measured values (Fig. 6a), except south of 55°S where the model values are too high in August. This is produced by the August λ(x) function (Fig. 5a), which has a maximum southward going survival distance at 60°S and a sharp reduction south of 60°. The vapour flux data net is sparse for the high southern latitudes and consequently the resulting λ(x) function could be in error south of 55°S. Also there is considerable uncertainty in the monthly evaporation of E fields poleward of 40° latitude. Poleward of 40° annual values of E are used. Figure 7 shows the August measured and modelled δs and ds and the model values do in fact deviate significantly from the data for some of the higher southern latitudes. Figure 8 gives the February results and data. There are clear seasonal variations in δ and d and with the above exception the model simulates their zonal averages. This is further demonstrated by Figure 9 which plots (δaug - δfeb) for the model and data versus latitude.
Fig. 6. Zonal averages of February (solid) and August (dashed) precipitation rates, (a) measured and (b) model.
Fig. 7. August values for: (a) δ(18O), (b) d. Dots indicate data, lines the model.
Regional Solutions
The regional solutions start at a given coast or sea-ice front (x = xc) using as input the λ(xc) and δs provided by the global sea-level model. The regional part of the model then calculates the δs and ds at target sites x t along some trajectory x t > x c over which precipitation P(x), air temperature T(x), and surface elevation Z(x) can be specified. The regional solution calculates the survival distance function from P(x) and P(xc) using Equation (8). P(xc) is the precipitation rate off the coast. The trajectory should ideally be perpendicular to the precipitation isopleths and along predominantly stable vapour-flux paths. If Z(x) is a constant close to zero (e.g. over an ice shelf) the sea-level integration procedure by itself would give reasonable δs. If, however, the trajectory goes over mountains and/or ice caps, it is necessary to also “model” the fact that for a given target site xt > xc and Z(xt) > 0 the more locally derived water vapour has not been mixed vertically as much as more distant source’s vapour.
For a given target x t beyond the coast with Z(xt) > 0 the regional solution starts with the sea-level value for the vapour contributed by the various sources at x, i.e. y(x,xt) given by Equation (1). The sea-level precipitable water
Fig. 8. February values for: (a) δ(18O), (b) d. Dots indicate data, lines the model.
Fig. 9. Seasonal δ(18O) difference. (δaug - δfeb) versus latitude. Dots indicate data differences, lines are from the August and February model runs.
W(x t) is given by Equation (2). One wants to reduce y(x,xt) for sources (xs), that are closer to xt. One also wants to preserve the known vertical structure for precipitable water W(xt,Z) (e.g. Eriksson, Reference Eriksson1965) which can be closely approximated by
where 2.0 < Lz < 3.0 km and where W(xt) is the sea-level value. This form of W(xt,Z) can be preserved and the more local water vapour excluded from higher-elevation target sites if in Equations (1) and (2) y(x,xt) is multiplied by exp{−(Z(xt)/Lz)Lx/(xt – x)} where Lx is the horizontal mixing length, i.e. that meridional horizontal distance needed to get a given source’s moisture mixed up to an elevation Lz. Thus for regional solutions Equation (2) is replaced by
The meridional horizontal mixing length Lx is estimated from the diffusive water-vapour flux terms TE + SE (Eriksson, Reference Eriksson1965) using [Inline 9] where K* is the meridional diffusivity (Austausch coefficient) for water vapour. Using 850 mbar values of TE and SE and q (Oort, Reference Oort1983), K* can be calculated as a function of latitude. It is not constant but varies in the range 0.18 × 107 < K* < 0.4 × 102 m2 s−1. The average mixing length is then estimated by [Inline 10] where [Inline 11] is the average meridional wind speed. In this work Lx = 2700 km in the Northern Hemisphere and 2000 km in the Southern Hemisphere. Equation (13) is not derived but it does accurately (to within 10%) reproduce the vertical precipitable water function and excludes the nearer sources in a plausible way. As is shown below, the results using Equation (13) are in agreement with other work.
Temperature inversion and precipitation weighting
For model runs over high-latitude non-marine regions it is often the case that the surface temperatures are lower than the condensation temperatures, i.e. there is a temperature inversion. This is particularly important in cold areas where the model must choose temperatures to feed to the Si(T) function if Τ < Ts. The Antarctic inversion function of Jouzel and Merlivat (1984) is used here.
For annual average runs the annual average temperature over cold non-marine areas is not likely to be equal to the year’s precipitation-weighted temperature. Usually the latter is several degrees warmer than the annual average. Annual regional solutions should use the precipitation-weighted temperatures. These are obtained by “correcting” annual averages using nearby meteorological station data.
An example of a regional solution
The snow accumulation map (from August to June) for the Queen Elizabeth Islands, Canada, is shown in Figure 10 (Koerner, Reference Koerner1979) along with a trajectory for a regional solution. The path starts at the coast xc and the Ρ function along it comes from the map values. The path’s elevation function Z(x) is taken from topographic maps up to Ζ = 1800 ma.s.l. after which Ζ is held constant. This is done because Koerner (1979) has shown that isotopie depletion beyond about 70 km from the coast is isobaric and that there is for (X - Xc) > 70 virtually no elevation effect in stable isotopes. The annual air temperature at 1800 ma.s.l. near the tops of the local ice caps is about −24° (Paterson and others, Reference Paterson1977; Fisher and others, Reference Fisher, Koerner, Paterson, Dansgaard and Reeh1983) and this is taken as a constant for (X - Xc) > 70 km. Up to 70 km, T(x) is found assuming a lapse rate of −0.61°/100m (Wilson, Reference Wilson1967). The global model provides the coastal start values for the δs and λ. The offshore precipitation P(xc) that is used as input to Equation (8) is obtained from Jaeger (1983) and local coastal station seasonal precipitation distributions, i.e. P(xt) = 9.3 g (cm−2 a−1). At the same time as Koerner mapped snow accumulation he also sampled that year’s snow for δ(18O). Measured δs versus distance from the coast along the trajectory in Figure 10 are shown in Figure lla as dots. The solid line in Figure 11 gives the modelled 8s assuming a temperature inversion and a precipitation weighting correction to annual air temperatures of +4.5°.
Fig. 10. Precipitation (cm of w.e.) from August 1973 to June 1974 in the Canadian Arctic islands and an example vapour trajectory perpendicular to the isopleths. The regional solution described in the text is along this trajectory. Map adapted from Koerner (1979).
The dashed line is a model run with no inversion and no precipitation weighting on temperatures. Figure 11b shows d the deuterium excess for both cases. As yet there are no d data for comparison. It is fortuitous that the global model provided a coastal δ very close to the actual measured one (Fig. 11a). In general the regional solutions could be expected to have an offset from the measured values. Figure 11a demonstrates that the model runs agree well with the data showing a steeper “orographic” depletion from the coast to (X - Xc) = 70 km and a gentler near-linear “distance from coast” depletion that is in good agreement with Koerner’s (1979) result.
Similar calculations have been made in regions where there is precipitation and δ data: north-west Greenland,
Fig. 11. Regional annual average solutions along the trajectory in Figure 10 versus distance from the coast, (a) Dots are measured δs (Koerner, Reference Koerner1979). The solid line is from the model assuming there is a temperature inversion over the islands and a correction to input temperatures due to seasonal weighting by precipitation rate. The dashed line is a model result with neither of the above assumptions, (b) Model deuterium excess d results for the two cases above. There are virtually no d data available for comparison.
Ross Ice Shelf, East Antarctica, Victoria Land, and mid-and south Greenland. In all cases except the last two in Greenland the model simulated the stable isotopes rather well from sealevel up to the highest elevations. These detailed local solutions will appear elsewhere. The mid- and south Greenland areas seem to be examples of areas that receive significant precipitation from along several vapour trajectories. In such cases it is necessary to add precipitation from several paths in order to model the measured δ pattern.
Differences in source zones for ice-cap sites in the two hemispheres
For a given target site at sea-level, or at any given elevation inland, the model calculates the per cent contribution of each source zone to the site’s precipitation. Figure 12 shows the source contributions per cent for 5° zonal strips for Crête (Greenland), a site at Crête’s latitude on the Greenland coast and for Vostok (Antarctica), this model predicts that the Crête site (at Ζ = 3100 m a.s.l.)
Fig. 12. The percentage contribution of 5° zones to precipitation at various sites in two hemispheres, (i) Crête, a site on Greenland at 3100ma.s.l. (ii) A coastal Greenland site at the same latitude as Crête, iii. Vostok, a station at 3488 m a.s.l. in East Antarctica. Antarctic sites evidently receive a higher proportion of precipitation originating from high-latitude cold oceans than sites in the more “continental” northern hemisphere.
receives over 70% of its accumulation from between 20° and 55°N. The predicted contribution weighted latitude for Crête is 35°N. Using a single source isotope model, Johnsen and others (1989) estimated that precipitation around Crête originates from oceans that have a surface temperature of about 20–25°C. This corresponds to a zone between 35° and 21°N.
The Vostok Station (Z = 3488 m a.s.l.) moisture seems to come from higher latitudes with maximum contributions from about 55°S. The predicted contribution-weighted latitude for Vostok is 44°S. Jouzal reports (personal communication) that vapour-tracing experiments with a GCM (including stable-isotope fractionation) gives the predominant source ocean temperature as 10°C for high East Antarctic sites. This temperature corresponds to latitude 43°S.
Thus, both this model and other work suggest that there is a significant difference in the latitudes contributing to high-elevation polar sites in the two hemispheres. This difference is reasonable because the Southern Hemisphere is mostly ocean between 30° and 65°S and the wind speeds are high. By contrast, the Northern Hemisphere is much more continental in these latitude zones and the wind speeds are lower (see Fig. If and Id). The large proportion of source moisture from cold southern oceans explains the deep coastal d minimum at 70°S (see Fig. 3b).
A model run for 75°N at sealevel and one at 1800 m a.s.l. (for the Devon Island Ice Cap. Arctic Canada) predict that 21% of precipitation at sea level originates north of 60°N and that at 1800 m a.s.l. the same source zone only contributes about 9%. Koerner and Russell (1979) on the basis of their isotopie measurements concluded that at sea-level on Devon Island 18–25% of the precipitation was locally derived; and that at the top of the ice cap only 8% came from Baffin Bay.
It would seem that Equation (13) and the regional model solutions give realistic results and allow the source zones to be identified.
Conclusions and Further Work
Using evaporation rates, precipitation, and vapour flux as the fundamental inputs to modelling the global and regional stable-isotope patterns seems to provide a relatively simple conceptual framework for explaining most of the observed features on both scales for annual average and seasonal data. In that precipitation rates in cold high-elevation areas are closely linked to vapour-mixing ratios and hence to air temperatures, the δ-air temperature relationship (Dansgaard and others, Reference Dansgaard, Johnsen, Clausen and Gundestrup1973) should be expected to hold as long as the distribution of source zones and the vapour-flux regimes remain constant. However, if there are major changes in global/local circulation patterns the present relation might break down. One object of future work will be to use existing long time series of stable-isotope records at various latitudes to try and extract the model’s input variables (i.e. run the model backwards).
The model can also in principle be used to calculate marine-derived impurities in precipitation and thus can be utilized to try and reproduce the interrelationships between the measured stable isotopes and these impurities.
