1. Introduction

The present paper concerns the trapping of air bubbles as firn (compacted but permeable accumulated snow grains) becomes non-permeable ice as the density increases. The bubble trapping is analysed using the percolation model from lattice statistics. The statistics of this air-bubble trapping are of some importance because the composition of air from bubbles in polar ice provides a record of atmospheric concentrations of climatically active atmospheric constituents such as CO₂. These records span periods ranging from decades to millennia, depending on the accumulation rate and thickness of the ice cap from which the record is obtained (see, for example, Reference Delmas, Ascencio and LegrandDelmas and others, 1980; Reference Neftel, Moor, Oeschger and StaufferNeftel and others, 1985; Reference Pearman, Ethridge, de Silva and FraserPearman and others, 1986; Reference Barnola, Raynaud, Korotkevich and LoriusBarnola and others, 1987). As well as being required for the direct interpretation of concentrations of trapped gases, the statistics of the firn-closure process are also important in the interpretation of those measurements of firn closure that are to be used in the interpretation of gas data.

The percolation model is a statistical model describing random networks. It has been used to model a wide range of physical processes (Reference Zallen, Deutscher, Zallen and AdlerZallen, 1983). The relevance of the percolation model to bubble trapping in firn was recognized by Reference EntingEnting (1985b) and Reference Stauffer, Schwander and OeschgerStauffer and others (1985). When modelling firn closure in terms of percolation, we can draw on a very large body of existing work on this and related models. Section 2 below reviews some key results from this field of “lattice statistics” and introduces some of the technical terms involved (see also the review by Reference EntingEnting (1987)). The present study is motivated by the fact that many of these lattice-statistics models, including the percolation model, exhibit transitions characterized by fluctuations on all length scales. Reference EntingEnting (1985b) identified the fluctuations in the trapped bubble volume, as measured by Reference Schwander and StaufferSchwander and Stauffer (1984), with the statistical fluctuations of the percolation transition. In order to investigate the significance of these fluctuations for gas dating, we performed Monte Carlo simulations using the ultramarine lattice model introduced by Reference Stauffer, Schwander and OeschgerStauffer and others (1985).

The percolation model, as used here, is a descriptive statistical model of firn closure, rather than being fully mechanistic. In particular, the model is expressed in terms of the probabilities of channel closure and does not attempt to predict these from the physical properties of the firn. This approach is motivated by the recognition that the statistical aspects are the most important characteristics of the firn-closure process. Any more comprehensive model that seeks to predict properties such as the close-off density must take the percolation statistics into account as a prerequisite for defining the close-off point. Some of the issues involved in modelling firn closure have been discussed by Reference StaufferStauffer (1982). The universality theory property of lattice model transitions (see section 2 below) predicts that, suitably scaled, a simple percolation model will correctly represent the characteristics of more complicated percolating systems.

One important qualitative result from this identification was the recognition that, because of the sharp cutoff in the trapping, the deconvolution of atmospheric composition from a record of mean bubble composition is a relatively well-conditioned inversion problem (Reference EntingEnting, 1985a). Reference Enting and MansbridgeEnting and Mansbridge (1985) (henceforth referred to as EM) suggested that, because of this relatively well-conditioned nature of the process of deconvolving bubble compositions to obtain a record of atmospheric composition, the trapping could be parameterized by an effective age. This effective age differed from the age of the ice by a delay time which was expressed as an integral over the trapping distribution (see Equation (6) below). From the simulations described below, we can estimate the sample-to-sample variation in this delay time.

The lay-out of the remainder of the paper is as follows. Section 2 reviews some key concepts from lattice statistics. Section 3 reviews the concepts of the bubble-trapping distribution and the effective age of the gas. A number of relevant measurements of ice-core properties are noted. Section 4 describes the lattice model used in the Monte Carlo studies. Section 5 presents the results of the simulations. Section 6 gives a brief discussion of the effect of pore closure on the restriction of diffusive mixing through the firn, and the role of percolation theory in understanding this effect. Section 7 gives some concluding remarks while additional details of the simulation algorithm are given in the Appendix.

2. Critical Phenomena and Percolation

In this section, we review some of the relevant results from lattice-statistics models of critical phenomena, with particular reference to the percolation model. Further details may be found in Reference EntingEnting (1987) and references therein. Some technical terms of particular relevance to the firn-closure problem are noted below in italics.

The term lattice statistics is used to cover a range of mathematical or statistical models defined at discrete points in space, usually on a regular lattice. The properties of the models are derived from locally defined statistical characteristics. Many of the applications have been in solid-state physics where the lattice corresponds to an actual crystal lattice and the statistical variations arise from the effects of temperature. These models have been the subject of very extensive study. The references in the review by Reference EntingEnting (1987) can be used as an initial point of contact with the lattice-statistics literature. In addition, Reference EfrosEfros (1986) and Reference StaufferStauffer (1985) have published introductory books specifically on the percolation model. The first of these covers a greater range of topics but suffers from the lack of any references.

The percolation model has a very simple statistical characterization. In the form used to model bubble trapping, the bubbles are associated with points on the lattice, the sites. Neighbouring sites are regarded as being linked by bonds on the lattice. These bonds are regarded as being either “connected” or “broken”. The connectivity is assigned randomly, independently for each bond. The model deals with the statistics of clusters of sites joined by chains of “connected” bonds. When modelling firn closure, the “connected” bonds represent open pores or channels between the snow grains.

There is a wide class of lattice models (including percolation models) which exhibit what are known as critical phenomena. The characteristics of these critical phenomena are:

• i. There is a critical probability at which the macroscopic properties of the system change abruptly. In the percolation model, there is a particular (lattice-dependent) probability, p _c, above which arbitrarily large connected clusters will occur, but below which the mean cluster size remains finite. In other lattice models, the abrupt change is used to model a change in thermodynamic phase as a function of temperature. The general term critical point is used to define the point at which the change occurs.

• ii. At this critical point, there is no characteristic length scale; the spatial characteristics show statistical self-similarity (i.e. various statistical properties remain invariant under a change of length scale). In the terminology of Reference MandelbrotMandelbrot (1977), clustering properties are described by fractals. Correlations between sites show power-law dependence on distance at the critical point.

• iii. Various other statistical properties exhibit singularities at the critical point. The behaviour in approaching the critical point can often be represented by power laws. The powers involved are called the critical exponents with standard symbols (usually Greek letters) denoting many of these exponents. For example, as a critical point is approached, characteristic length scales (denoted ξ) diverge with a power-law dependence on the distance from the critical point: Such power-law scaling relations apply within a critical region surrounding the critical point.

• iv. The values of the critical exponents are independent of many local details of the statistical model. In particular, for percolation models, the critical exponents are the same, whether one considers clusters of randomly occupied sites, clusters of randomly occupied bonds or, as in the model described below, clusters of sites connected by randomly occupied bonds. The critical exponents are also independent of the particular lattice, apart from depending on the dimensionality. These independence properties are referred to as universality. An extreme example is the use of a lattice model to represent the liquid-gas co-existence line and the liquid-gas critical point that terminates this co-existence line. This model can represent many features of such critical points even though neither phase occurs on a lattice.

• v. As the critical point is approached, various statistical quantities become extremely sensitive to small perturbations. In the liquid-gas model, the compressibility becomes arbitrarily large. In magnetic systems, the susceptibility (the rate of change in magnetization as the external field is increased) diverges. In the percolation model, there is an analogous divergence in the mean number of sites per connected cluster.

• vi. As well as the large responses to external forcing, the statistical properties near the critical point exhibit large intrinsic fluctuations. In the liquid-gas system, these correspond to density fluctuations that are observed through anomalously high light scattering, known as critical opalescence. Similar statistical fluctuations occur in the percolation model. Reference EntingEnting (1985b) interpreted the highly variable results for the amount of trapped gas near the firn close-off (as measured by Reference Schwander and StaufferSchwander and Stauffer (1984)) as an example of such critical fluctuations. These observations were a key piece of evidence in arguing for the applicability of the percolation model in representing bubble trapping in firn. These critical fluctuations are an intrinsic part of the transition. Indeed, critical phenomena are sometimes referred to as co-operative phenomena because the transitions occur due to the cooperative effects of such fluctuations.

The most common forms of percolation model are site percolation (clusters of randomly occupied sites where pairs of occupied neighbours are regarded as connected) and bond percolation (clusters of randomly occupied bonds). Reference StaufferStauffer and others (1985) modelled firn closure using a conventional bond percolation model. Reference EntingEnting (1986, Reference Enting1987) proposed that, to represent bubble trapping, a more appropriate form of the percolation model would consider clusters of sites connected by randomly occupied bonds. In this “bubble-trapping” percolation model, the sites represent space between ice grains and the bonds represent connecting channels that are gradually closed off. A key statistic in percolation modelling is the percolation probability, P, which is the probability of a site or bond being part of an infinite connected cluster. At the critical point, Ρ goes to zero as with β ≈ 0.454 in three dimensions.

Reference EntingEnting (1986, Reference Enting1987) noted that the bubble-trapping model (clusters of sites connected by randomly occupied bonds) was the special form of a percolation model obtained by taking the q → 1 limit of the q-state Potts model (Reference Fortuin and KasteleynFortuin and Kastelyn, 1972; Reference WuWu, 1978). This connection showed the relation between the percolation model and more general statistical mechanics models of phase transitions. It was also useful in the derivation of-series expansions for the trapping distribution (Reference EntingEnting, 1986), but otherwise gives little insight into the percolation model and so will not be considered further in this paper.

Reference EntingEnting (1987) listed five main classes of mathematical technique that have been used to investigate critical phenomena in lattice models. These are (a) exact solutions (generally not available), (b) statistical closures (usually inaccurate in the critical region), (c) power-series expansions of statistical moments (as applied to firn closure by Reference EntingEnting (1986)), (d) renormalization group techniques and (e) Monte Carlo simulations.

Monte Carlo simulations are used in the present study because the statistical fluctuations are simulated directly. However, they are restricted by the fact that the simulation can only be for a finite system. The other computational techniques above allow a formal limit of arbitrarily large systems to be considered. Since the range of spatial correlation diverges near the critical point, finite-size effects become important and must be corrected for.

An example of the importance of correcting for finite-size effects is that Reference EntingEnting (1986) estimated the p_c = 0.39 ± 0.01 using series expansions based on combinatorial factors that involved less than 20 sites, but which were not subject to edge effects. In contrast, the Monte Carlo calculation by Reference StaufferStauffer and others (1985) estimated p_c = 0.42 on the basis of simulations involving over 4000 sites. The results are presented below, involving up to 165 888 sites, show that the apparent p_c decreases with sample size and appears to be converging towards a value in the range estimated from the series calculations.

Since the main object of the present paper is an analysis of the size-dependence of the sampling statistics of bubble trapping, it is necessary to distinguish this from the finite-size effects noted above. Finite-size effects inevitably occur in Monte Carlo simulations because only a finite number of sites can be simulated. Within such finite simulations, we consider the question of how the bubble-trapping statistics are affected by taking specific sized samples for gas analysis. In order to study this second type of finite-size effect, the statistics are taken from a sub-region of the sites simulated.

3. The Bubble-Trapping Distribution

In the analysis of bubble trapping, there are a number of alternatives for the independent variable. The percolation models work in terms of the probability of bond closure (u) while levels within an ice core can be described in terms of density ρ, depth (d) or age (z). In ideal situations, the four variables u, ρ, d and z are all monotone increasing functions of each other. Results reviewed by Reference Schwander and StaufferSchwander and Stauffer (1984) suggest that the critical probability u_c, beyond which trapping is complete, corresponds to ρ_c = 0.82 M gm⁻³ under a range of different deposition rates. There is, however, a dependence of closure density on the temperature of deposition (personal communication from D. Raynaud). Further studies of aspects of firn closure have been given by Reference StaufferStauffer (1982), Reference LoosliLoosli (1983), Reference Stauffer, Schwander and OeschgerStauffer and others (1985), Reference Schwander, Stauffer and SiggSchwander and others (1988) and Reference Schwander, Oeschger and LangwaySchwander (1988).

The relation between depth density and the age of the ice depends on the local snow-accumulation rate which varies by several orders of magnitude between different sites. The relationship between depth and age can be established in a variety of ways. In regions of high accumulation, seasonal variations in quantities such as electrical conductivity or oxygen-isotope ratios can be used to count the number of years since deposition. In regions of slower accumulation, ages must be inferred indirectly by using glaciological models.

EM (i.e. Reference Enting and MansbridgeEnting and Mansbridge, 1985) described the basic concepts involved in interpreting the mean concentrations of gases recovered from ice-core samples. They defined R (z,z_t ) as the amount of gas that was trapped z_t years ago in ice that was deposited z years ago. Clearly, this is zero for z_t > z, i.e. the gas cannot be trapped until the ice is present. The existence of a sharp transition between firn and ice implies that there is also a lower bound on z_t. The gas must be trapped within a time of z_c of the snowfall. Thus, for non-zero R, we require

(1)

If we assume a uniform deposition process, then we set

(2)

The sharpness of this trapping distribution determines the extent to which gas bubbles from different times are trapped in a single layer. The cumulative trapping distribution was defined by EM as

. (3)

They used this to define the normalization by requiring

. (4)

EM showed that, if ambient gas concentrations c(z) were changing linearly, then the average concentration, q(z), occurring in bubbles in an ice layer of age z corresponded to ambient concentrations at a younger age, z_eff, so that

. (5)

with z_eff defined in terms of the mean delay or mean trapping time z by:

(6)

In its simplest form, neglecting finite-size effects, the percolation model of bubble trapping equates 1 − S(z) with the percolation probability, P. For an infinite lattice, percolation theory defines Ρ as the probability that a given site is part of an infinite cluster of sites connected by open channels. For finite samples, the corresponding quantity will exhibit statistical fluctuations as shown in the simulations presented below. In view of the analysis of gas ages summarized above, the most important question concerning critical fluctuations is the size of the critical fluctuations in the mean trapping time z . It should be noted that since, by definition,

(7)

we have

(8)

independently of z_x.

When dealing with finite samples, the point of complete trapping will show fluctuations about the infinite size limit z_c. For this reason, z_x must be chosen to be greater than the actual point of complete trapping for any sample under consideration.

The analysis summarized above describes concentrations, q(z), for layers in which the air is completely trapped. For concentrations, q(z), obtained in ice layers that are not completely closed off from the atmosphere, the expressions above must be modified. What is required (see EM) is to rescale S(z) so that it describes incomplete trapping. The age, z, of the ice layer replaces z_c in Equation (6) to give

. (9)

This aspect of the model will not be considered in the present study.

4. The Lattice Model

The application of the percolation model to firn closure is based on the following propositions:

• i. The snow grains of the firn form a random network of spaces and channels containing air.

• ii. As the firn is compressed, the narrowest gaps are closed, leaving bubbles of air in the larger spaces. A more detailed analysis of the changes occurring during firnification has been given by Reference StaufferStauffer (1982).

• iii. The closure of channels will be effectively a random process reflecting the random structure of the firn.

• iv. The percolation model is the appropriate tool for analysing the statistics of random closures. The various properties described in section 2 above will apply.

• v. Percolation modelling predicts a sharp transition between the state where the fim is permeable and the impermeable ice. This transition reflects the long-range statistics of connectivity and will occur before all channels are closed.

• vi. The universality property of percolation models predicts that the main features of the bubble trapping at the firn-to-ice transition will not depend on the details of the model. The only efTect of changes in such details will be a smooth rescaling of model parameters. Among such “irrelevant” details are: (a) the use of a regular lattice rather than a random network, (b) local correlations in channel closure, (c) assigning all of the trapped volume to the larger spaces that become bubbles and ignoring the air in the channels.

• vii. While universality enables us to relate critical behaviour on different lattices, the comparisons will be most direct if the choice of lattice reflects the statistics of connectivity in the firn. On the basis of measurements of firn structure (see Reference StaufferStauffer, 1982), Reference Stauffer, Schwander and OeschgerStauffer and others (1985) suggested that an appropriate lattice model for bubble trapping could be obtained by taking the ice grains as being on a body-centred-cubic (b.c.c.) lattice.

• viii. The statistics will exhibit large fluctuations near the firn-to-ice transition.

• ix. If the trapping times are sufficiently long compared to the times of mixing through the firn, then the age distribution of the trapped gas will correspond to the age distribution of the time of trapping.

• x. The percolation model considers the percolation probability, P. This is the probability that a given site will be part of an infinite cluster of connected sites rather than part of a finite cluster. The bubble-trapping model regards sites in finite clusters as being cut off from the surface and so the proportion of trapped bubbles is 1 − P.

• xi. The function 1 − P(u) gives the mean proportion of trapped bubbles in an infinite system. Reference EntingEnting (1986) applied series-expansion techniques in estimating this function. The present paper considers the extent to which the proportion of trapped bubbles in finite samples will show a sample-to-sample variability about the mean.

• xii. The most appropriate way to study the statistics of bubble trapping is through Monte Carlo simulation of a limited region of firn.

Following the earlier work by Reference StaufferStauffer and others (1985) and Reference EntingEnting (1986), we model the firn by locating snow grains on a b.c.c. lattice with the interstitial sites representing the locations of air bubbles. The sites each are joined to their four nearest neighbours by bonds that are connected (with probability p) or broken (with probability u = 1 − p). We describe the lattice in terms of cubic cells containing 12 sites and 24 bonds as shown in Figure 1, and consider arrays of m × m × n of these cubic cells. We consider the statistics of connection to the upper boundary (layer 1) as an approximate representation of long-range connection to the surface. We apply periodic boundary conditions in the horizontal (m × m) layers, i.e. each horizontal cross-section is taken as being surrounded by identical copies of itself. The simulation procedure assigned each bond (connecting sites i and j) a randomly chosen “cut-off time”, u_ij which corresponded to the “time” at which bond closed. The u_ij are described as “times” in order to give a more intuitive description but are more correctly regarded as the probabilities of site trapping and bond closure and therefore the u_ij were selected uniformly from the range [0,1].

Fig. 1. The basic simple cubic cell and the 12 sites and 24 bonds defining the lattice of b.c.c. interstitial sites. Each site lies on one of the square faces defined by two of the axes and two of the dashed lines. Dashed sites are regarded as lying in cells neighbouring the one shown.

The model statistics of bubble trapping were obtained from the site-closure values u_i, the “time” beyond which site i was no longer connected to the surface, i.e. the “time” at which the last path of open bonds connecting site i to the surface was closed at some point. These u_i were calculated from the u_ij using the algorithm described in the Appendix. The bubble-trapping statistics were computed using a sub-set, A, of sites so as to reduce the influence of the finite size of the simulation. Specifically, the bubble-trapping distribution was

(10a)

where

. (10b)

Consequently X(u) is the number of sites cut off from the upper boundary (i.e. trapped) at “time” u and S(u) is the proportion of trapped sites at “time” u.

The lattice model involves several key approximations:

• i. It regards all the trapped gas as being in the sites rather than being partly in the bonds.

• ii. It assumes independent closure probabilities for each bond.

• iii. It assumes a uniform lattice rather than a random structure.

• iv. It applies only to uniform deposition conditions.

The most important of these approximations is the assumption of uniform snow-deposition conditions. This is because, as noted in section 2 (point v), co-operative phase transitions such as the percolation transition, are characterized by divergent susceptibilities to external perturbations as well as divergences in the intrinsic fluctuations. The analysis of externally forced variations is largely outside the scope of this paper, except to note that changing conditions will mainly act to change z_c by changing the relation between z and u. Since the infinite size limit of S(z) has an infinite gradient at z_c, point-wise values of the trapping distribution near z_c can have finite responses to arbitrarily small perturbations. However, a change, Δ(z_c), in the critical point will induce only a change of the same order in the integral (8) defining the mean delay z .

The first three model approximations listed above will limit the quantitative accuracy but not the qualitative results of the present study. Universality arguments indicate that the critical exponents, such as v (which describes correlation lengths and finite-size effects) and β (which describes P(u)) will not be affected by such local details of the model. The amplitudes involved in the various critical phenomena will depend on the structure of the lattice. However, most of this dependence will involve gross features such as the mean number of neighbours for which the present model seeks to match observational data as closely as possible.

For the purposes of the simulation, the u_i and u_ij were discretized at intervals of 0.0001 and represented by 16 bit integers. The plots of S(u) are based on values tabulated as intervals of 0.001. In relating our simulations to actual ice cores, we have used estimates of 260 bubbles per cubic centimetre based on analysis of photographs of sections from the “ice” region of the BHD core analysed by Reference Schwander and StaufferPearman and others (1986). We have analysed our results as though the closure probability u was directly proportional to z. This is interpreting the intuitive description, given above, of the u_i and u_ij as times as literally true but it is a fairly crude approximation. The series-based estimate of the percolation probability Ρ on the lattice of b.c.c. interstitial sites was obtained by Reference EntingEnting (1986, Reference Entingh1987). The estimate 1 − Ρ is plotted as the solid curve in Figure 2 and suggests that the change from 10% closed off to 90% closed off occurs over a period of while for the BHD core the change from 10% to 90% closure occurs over about Δz ≈ 14 years ≈ 0.2z_c.

Fig. 2. Realizations of S(u), the proportion of trapped bubbles, for finite samples of bubbles, plotted as functions of bond-closure probability, u. The smooth curve in each case is the infinite lattice estimate obtained from series expansions. Three or four realizations are shown in each case; the small roman numerals are a guide to distinguishing the lines. Some cases are only shown in part because of a high degree of overlap with other cases.

(a) Case 1: sample of 324 sites from lattice of 4752 sites; (b) Case 3: sample of 20 736 from 290 304 sites; (c) Case 5: sample of 324 sites from 36 228 sites.

In order to relate the simulations more closely to actual ice samples, the points of 10% and 90% closure give a better basis for matching the scales of age z and closure probability, u. Thus, ranges expressed in terms of u should be converted to time ranges by multiplying by Δ z/Δu with Δz determined for the particular core.

We have performed a sequence of simulations for various sizes of lattice and the sample sub-set A. As noted in section 2, there are two different types of finite-size effect involved. The first of these, which has been most extensively studied by other workers, is the effect of using a finite size of lattice for the simulations. For the present study, this effect represents a source of error that must be taken into account when analysing the simulations. The second finite-size effect, in which we are primarily interested, is the degree of statistical variability that arises from using finite samples of bubbles. This is a question of sampling statistics, a topic that has not often been considered in lattice statistics. In order to estimate the importance of the two finite-size effects separately, we have varied the lattice size and the sample size independently, considering the various cases listed in Table 1.

Table 1. Details of sizes of lattices and samples used in the simulations. The lattices always had a square cross-section with periodic boundary conditions in the horizontal plane. The samples used to determine the statistics were always cubes. All dimensions are given in terms of the cubic cell containing 12 sites, except that the sample size is given as the number of sites

The samples have their linear dimensions increasing by a factor of 2, and involve 324, 2592, 20 736 and 165 888 sites. For comparison, we estimate that the 1 kg samples analysed in our laboratory contained about 290 000 bubbles (on the basis of 260 bubbles cm⁻³ and a density of 0.9 Mgm⁻³). Reference Stauffer, Schwander and OeschgerStauffer and others (1985) estimated that the 35 cm³ samples used in their work and that of Reference Schwander and StaufferSchwander and Stauffer (1984) contained about 30 000 bonds (i.e. 15 000 sites on the ultramarine lattice). Earlier studies by the Bern group used samples as small as 1 g which would correspond to about 400–500 bubbles.

5. Results

The three parts of Figure 2 each show four realizations of S(z) for cases 1, 3 and 5 of Table 1. Also shown, as the smooth line on each figure, is the estimate of S(z) derived from the series expansions given by Reference EntingEnting (1986). The estimate is expressed as 1 − [V₈(u)/W₇(u)]^0.454 where V₈(u) and W₇(u) are polynomals of degree 8 and 7, respectively. Specifically, V₈(u)/W₇(u) is the Fade approximant to P(u)^1/0.454 obtained from the series for P(u) given by Reference EntingEnting (1986).

Figure 2a and b (cases 1 and 3) illustrate a systematic tendency (confirmed by the results of cases 2 and 4) for the smaller samples to produce small values for the cut-off probability so that, as the sample-size increases, the cutoff tends towards the series estimate. This size dependence is primarily due to the finite size of the lattice used. Comparing Figure 2a and c shows that the cut-off probabilities for the smallest samples tend to increase as the size of the surrounding lattice increases. The degree of variability between the realizations remains essentially unchanged. This size dependence of the average cut-off seems to explain much of the discrepancy between the values of u_c, estimated from series techniques as 0.61 ± 0.01 (Reference EntingEnting, 1986), and from Monte Carlo simulations of up to 8748 bonds, as 0.58 (Reference Stauffer, Schwander and OeschgerStauffer and others, 1985).

The various parts of Figure 2 show the variability (and size dependence of this variability) in the cumulative trapping distribution as a whole. As noted above, in studies of gas trapping, we are primarily concerned with the mean trapping time, z defined by Equation (6), and wish to determine how much sample-to-sample variation can occur in this quantity. It is to be expected that the mean trapping time will be less variable than the various aspects of S(z) because it can be defined as an integral over S(z).

However, it must be emphasized that for any single sample, the deviations between S(z) and its ensemble mean are not independent for various z values. In particular, the deviations must be such that S(z) → 1 for sufficiently large z. Indeed, the need to construct quantities based on the closure history of single samples was an important reason for using Monte Carlo simulation in this study. Alternative approaches such as using statistics of cluster size to determine finite-size effects would produce ensemble averages at each u value rather than statistics describing the closure histories of individual samples.

Figure 3 shows the cumulative distributions of the mean trapping times based on 20 realizations for the three smallest cases. It will be seen that, except for the smallest (324 site) examples, the distribution of trapping times is small compared to the mean time over which closure takes place. These results indicate that the mean closure probability is relatively insensitive to the critical fluctuations that affect finite-size samples.

Fig. 3. Cumulative distributions of mean trapping times obtained from 20 realizations each of cases 1, 2 and 3. Solid curve, case 1; dashed curve, case 2; dotted curve, case 3.

Finally, we consider the scaling properties of the distribution of trapping times shown in Figure 3. For 20 L × L × L samples, we order the 20 mean trapping times z _j and, as a robust measure of the spread, consider the interquartile range . These are shown on a logarithmic plot against L in Figure 4. Finite-size scaling would suggest τ∼L ^–1/v with v ≈ 0.9 in three dimensions. The scaling line is shown on Figure 4 and is consistent with the data, although our samples are too small to give a precise estimate of v What is clear is that T is decreasing far more slowly than the L^−D/2 (with D = 3) that would occur for independently trapped bubbles. This emphasizes the co-operative nature of the bubble-trapping process.

Fig. 4. Interquartile range of the mean trapping time plotted against sample length L. The results were from the sets of 20 cases used to produce Figure 3. The scaling line aL^−t/v with v = 0.9 is shown dashed.

6. Extensions

The analysis above has described the effect of finite sample size on the statistics of bubble-trapping time. Transferring these results to the statistics of gas age is only valid when the trapping time is long compared to the mixing time through the firn. For cores from areas of even moderately high snow accumulation this condition will cease to be valid.

Reference Schwander, Oeschger and LangwaySchwander (1988) has performed some modelling studies of the cases where the effective gas age is influenced by the slowness of gas mixing through the firn. The firn properties were determined by laboratory measurements. The percolation model can shed further light on this problem of diffusive mixing through the firn near close-off.

The relevant problem in percolation theory has been studied in terms of the problem of conductivity of random resistor networks. The conductivity is the analogue of the diffusivity in a network of channels. The conductivity problem in percolation theory is a matter of some difficulty. There is apparently no series-expansion formalism for calculating the conductivity, although related statistical moments can be calculated and the critical behaviour of the conductivity can be inferred from these. It is also possible to determine the conductivity from Monte Carlo experiments, either by the tedious “brute-force” method of solving Kirchoff’s equations on a set of random networks, or more elegantly through a “self-averaging” tiansfer-matrix method described by Reference Derrida, Zabolitzky, Vannimenus and StaufferDerrida and others (1984). The result is that, as expected, the conductivity goes to zero at the critical point. Most importantly for the ice-core application, the critical exponent is greater than 1 so that the gradient of the conductivity curve is zero at the critical point, in contrast to S(z) where the gradient is infinite. This means that the conductivity in a resistor network (and diffusivity in a network of channels) may be small throughout a significant fraction of the critical region and possibly beyond.

These properties of the conductivity/diffusivity problem have important implications for gas dating, depend-ing on the relative time-scales of diffusion and bubble closure. If the bubble trapping is slow relative to diffusive mixing, then only the trapping distribution need be considered. For somewhat faster trapping, the effect of restrictions on the diffusive mixing will be to produce an age distribution that is more rounded than the trapping distribution. This will reduce the extent to which changes in atmospheric composition can be resolved from the ice-core record.

The most significant effects of restrictions on diffusion occur when differences between summer and winter accumulation lead to “winter” layers sealing off the “summer” layers below them even though the “summer” layers may be still permeable and have only a small fraction of bubbles trapped by local closure. In such a case, it would be expected that the statistical variability of trapping in “summer” layers is largely due to variations in the overlying winter layer and so much of the variability will involve year-to-year differences rather than within-year differences. The sample size becomes a relatively unimportant issue. In contrast, within the “winter” layers where the effective closure is due to changes within that layer, local statistical fluctuations are likely to play a much more important role.

7. Conclusions

The simulations described in this paper were carried out in order to assess the influence of critical fluctuations in the firn-closure process on the dating of gas samples trapped in polar ice. Since it has been possible to perform simulations of systems whose size is comparable to those of actual interest, it has not been necessary to use finite-size scaling to extrapolate the results of the simulation.

The main result of these studies is that the mean trapping time z is relatively independent of sample size over the range of interest. Furthermore, over this range, the sample-to-sample variation in the z is largely confined to a range that is small compared to the time over which most of the closure (firnification) takes place, except for small samples of order 1 g. However, the simulations also suggest that for gases with rapid changes in atmospheric concentrations (and for which the mean trapping time is not applicable) critical fluctuations may significantly affect the trapped gas concentrations.

The simulations suggest that the one important way in which critical fluctuations will affect gas dating is through the process of relating closure properties to actual ages. The conventional specification of this relation is to determine times at which 10% and 90% of the trapping has occurred. The simulation results shown in Figure 2 indicate that there can be a great deal of sample-to-sample variation in these two times, arising purely from the inherent statistical variability of the trapping process. The variation in the 90% trapping point is larger than the interquartile range for the respective sample sizes. Near the u values corresponding to 10% trapping, there is much less variability in the trapped volume but the flatness of the trapping distribution makes the position of the 10% point rather uncertain. This “calibration” problem will be exacerbated by the fact, noted above, that the critical region is one in which the system will be highly susceptible to the influence of external perturbations. As discussed above in connection with the effect of external perturbations, an uncertainty in the timing of the trapping due to the “calibration” will lead to pointwise values of S(z) being subject to large uncertainties. It turns out that z does not reflect this sensitivity; the uncertainty remains comparable to the uncertainty in the timing of the trapping.