Structural and Chemical Heterogeneity in Ancient Glass Probed using Gas Overcondensation, X-ray Tomography, and Solid-state NMR

Rare ancient glasses have complex, multi-scale structures requiring more sophisticated and non-destructive pore characterisation techniques than usual. Homotattic patch models for nitrogen adsorption gave better fits to the isotherm data, more accurate void space descriptors, and also greater understanding of the underlying physical factors affecting adsorption, than standard BET. These homotattic patch models revealed the critical role of iron impurities in determining adsorption behaviour. Non-destructive sodium-23 NMR relaxometry validated the homotattic patch model for some natron glasses, and, in turn, was validated using multiple quantum magic-angle spinning (MQMAS) Na NMR. X-ray tomography images of the glasses showed the presence of large macroporous bubbles, while FEG-SEM revealed nanopores within the glass matrix. A newly-developed, gas overcondensation technique, suitable for small amounts of low porosity material, assessed the inter-relationship between the disparate levels in this hierarchical porosity. This technique demonstrated that the nanoporosity did not form a ‘corona’ around the bubbles, due to leaching from the glass, as initially supposed from tomography data, but was completely disconnected, and, thus, is probably associated with glass alkalinity. Gas overcondensation is demonstrated as a non-destructive alternative to mercury porosimetry for probing multi-scale porosity in rare artefacts.


INTRODUCTION
Textural properties of amorphous materials, such as glasses, depend upon their fabrication method. These properties can thus be used as fingerprints to characterise particular types of materials, or the results of particular processing. For example, the surface area of amorphous materials may decline due to sintering occurring during calcination. It is thus necessary to obtain accurate measures of key textural properties using methods that can cope with the high level of chemical and structural complexity present in amorphous materials and still deliver accurate characteristics.
Specific surface area is a very commonly obtained descriptor of porous media, since it represents a key measure of the level of porosity present. The most common method used is the Brunauer-Emmett-Teller (BET) analysis of nitrogen adsorption isotherms [1]. Indeed, the International Standards Organisation make recommendations (BS ISO 9277:2010) concerning the experimental method and data analysis procedure to determine surface areas [2]. However, for complex, amorphous materials, like ancient glass, this standardised procedure may not be the most appropriate to deliver accurate results. The standard BET example, the polymer template in some ordered porous silicas gives rise to a microporous 'corona' around the mesopores [4]. Unfortunately, there are often limitations on the range of length-scales that can be probed together by standard methods. For example, imaging methods often have limits on the fineness of the pixel resolution and size of the field-of-view possible, whereas conventional gas adsorption experiments cannot probe macroporosity, and mercury intrusion porosimetry may be prevented by pore shielding due to narrow necks less than ~4 nm. However, gas overcondensation can provided statistically-representative information on macropores that are too large to detect using conventional gas adsorption experiments, but which are also shielded by pore necks too small to be penetrated by mercury intrusion [8,9]. This potentially removes any ambiguity about why macroporosity, known to exist within a given sample by other means, does not feature in either gas adsorption or mercury porosimetry pore size distributions. Further, gas overcondensation is not a potentially destructive method like mercury porosimetry (due to the high pressures involved), since the samples can be recovered undamaged or uncontaminated afterwards. This is important for rare or unique samples, such as archaeological finds. In this work nitrogen overcondensation will be used to probe the spatial disposition of nanoporosity relative to the very large macroporosity within the glasses. In this way, overcondensation experiments will test the hypothesis that nanoporosity was formed by leaching of species from the glass matrix surrounding the macroporosity.
In this work it will be shown that the standard BET analysis poorly accounts for the nitrogen sorption data and delivers inaccurate surface areas, and, in contrast, that the homotattic patch model more successfully accounts for adsorption on ancient glasses. It will be shown that the homotattic patch model alone reveals the importance to adsorption of the iron impurities in ancient glass. It will be seen that sodium-23 nuclear magnetic resonance relaxometry can independently validate this homotattic patch model for adsorption, and show that the model parameters have physical meaning. The presence of the different sodium environments determined from the T2 relaxometry will be validated by triple quantum magic-angle spinning solid state (3QMAS) NMR. Solid state magic angle spinning (MAS) NMR is routinely utilized to determine local environments in materials with limited long-range order, such as glass [10]. As sodium typically exists as a counter-ion in glassy materials then the range of different chemical environments, and hence the NMR chemical shift and quadrupole parameters, sodium experiences is small [11]. Therefore, significant resolution enhancement is needed to deconvolute the NMR resonances and to accurately determine the number of sodium environments in the glasses. Two-dimensional 3QMAS experiments separate the anisotropic quadrupole broadening (x-axis, F2) experienced by spin >1/2 nuclei, such as sodium, and the isotropic chemical shift information (y-axis, F1), allowing greater resolution to be achieved and overlapping sites to be deconvoluted [12]. However, the MAS technique is destructive as the sample needs to rotate at 12-13 kHz (720,000 to 780,000 rpm, to reduce anisotropic dipolar interactions between spins) and is not routinely applied to ancient artefacts.

Gas sorption
Common isotherm equations used to describe adsorption behaviour on porous solids include Henry's law [4]: where V is the amount adsorbed, Vm is the monolayer capacity, H is the Henry's constant, and P/P0 is the relative pressure. Henry's law is the low pressure limit for the Langmuir and BET equations mentioned below, and corresponds to the situation of sparse adsorption with no influence from adsorbate-adsorbate interactions, nor any indication of an ultimate adsorption capacity. Another common isotherm type is the Langmuir isotherm [4]: where B is the Langmuir constant describing the strength of interaction with the surface. The Langmuir equation is derived from a physical model of monolayer adsorption on a flat, homogeneous surface, but the mathematical form of the equation also corresponds to that found empirically for continuous micropore filling. In both the Henry's law and Langmuir equations the interaction constant is related to the heat of adsorption by the relation [4]: where E is the heat of adsorption, R is the gas constant, and T is absolute temperature. The pre-exponential factor is generally taken as approximately unity.
Real surfaces are often rough on the molecular scale. Concavities in the surface mean that the space for adsorbing the second and subsequent layers of an adsorption multi-layer declines with distance from the surface, and the maximum molecular capacity of each successive layer decreases. For surfaces that possess the special property of self-similarity over several lengthscales, and are, thus, known as fractals, the decline in the capacity of each layer of adsorbate is given by the equation [13]: where A1 is the area in the first adsorbed layer, Ai is the area in the ith adsorbed layer, and D is the surface fractal dimension (2D3). This effect can be incorporated into the BET model such that a fractal version of the BET equation is obtained. This is given by [6]: where C is the BET constant, and x is the relative pressure.
The effect of the fractal roughness, as described by Equation (4), is to lead to a decline in the amount adsorbed in each successive layer of adsorbate. Previous work has shown that a similar effect also arises when adsorption occurs as multiple, isolated, pyramidal molecular stacks, but the area for adsorption then obeys a more general power law such that [3]: where the power,  depends on the size of the base of the stack. For certain stack sizes the form of the upper part of the multilayer adsorption isotherm would, thus, be similar to that for adsorption on fractal surfaces. For example, adsorption as a pyramidal stack of base sidelength equal to 50 molecules on a flat underlying surface would appear the same as adsorption on a rough surface with fractal dimension of 2.24. Further, both effects can also occur together when adsorption occurs in stacks on a fractally-rough surface. Adsorption in isolated stacks is associated with patchwise chemically-heterogeneous surfaces.
The homotattic patch model was introduced in order to account for chemical heterogeneity [5]. This model considers the surface of the adsorbent to consist of a patchwork of different types of site each with their own characteristic adsorption behaviour. The model assumes that each of these patches is large, such that edge effects, where they neighbour other patches, are negligible. The resulting overall adsorption is thus a composite of the behaviour of the set of patches, such that: where Ii is the isotherm equation describing adsorption on the ith patch, and pi is the fraction of the surface occupied by patches of type Ii, such that the various pi-values obey: Capillary condensation in larger pores is predicted by the Kelvin equation [4]: where P/P0 is the relative pressure at which condensation occurs in a cylindrical pore of radius rp, k is a geometry parameter and depends on the pore type (for a cylindrical pore open at both ends k = 1; and for a pore with one dead end, or for desorption from a hemi-spherical meniscus, k = 2), γ is the surface tension and VM is the molar volume of the condensed liquid phase, θ is the contact angle with which the liquid meets the wall, and T is the absolute temperature. It is generally assumed, for nitrogen, that the adsorbed condensate is perfectly wetting of the surface such that the contact angle is zero, and thus the cosθ term is unity [4].
However, some studies [14] have allowed the cos θ term to be a free-fitting parameter, and, by doing so, it has been shown that a corrugated, cylindrical pore geometry can give rise to all of the IUPAC standard hysteresis loop types [15], thereby suggesting differences in porewetting can explain differences in isotherm shape. From Equation (9), the condensation and evaporation pressures for the same adsorbate in the same pore, or the condensation pressures for the same adsorbate in open/closed pores of the same diameter, or identical pores with different wetting properties, can be related via the ratio [16]: where the subscripts 1 and 2 refer to either condensation and evaporation, respectively, or two different pores of the same radius. For condensation and evaporation, for a through cylindrical pore with a fully wetting surface, k1 = 1, k2 = 2, and the cosθi terms both equal unity. Hence, in that case, the power  is equal to 2, and the relative pressure for evaporation is the square of the relative pressure for condensation. This case corresponds to the wellknown Cohan equations [17]. For less wetting surfaces, cosθi would be less unity, and, thus, the power would be less than two. Previous work, has shown that the power  to superpose adsorption and desorption branches for NLDFT kernels calibrated with fumed silica is 1.8 [16]. For a fully wetting equilibrium sorption system, with no hysteresis, the power would be unity.

SAXS
The intensity of X-ray radiation scattered by a fractal surface is known to be proportional to a negative power of the q vector, such that [18]: where q = 4πλ -1 sin(φ/2), λ is the wavelength of the radiation, and φ is the scattering angle.
Normally this behaviour is only observed if q satisfies the inequality q, where  is the characteristic length-scale for the structure creating the scattering. From the value of η, the fractal nature of the system under investigation can be determined [13]. If the exponent is in

Glass samples
The samples from Tell Brak, northern Syria (denoted BRAK 12, 14 and 17) date to the late 15 th and 14 th centuries BCE [19]. Their electron microprobe analyses (Table 1) show that they are soda-lime-silica glasses characterised by elevated levels of magnesium and potassium oxides when compared with the rest of the samples (below). These oxides are indicators that the soda source used to make the glass, probably in the form of sodium carbonate, was the ashes of halophytic plants. It is difficult to be more precise than this, but a possible genus used would have been Salsola [20]. The colorants used as indicated by the chemical analyses were a combination of copper to produce the colour and calcium antimonate crystals to produce opacity (BRAK 12), cobalt to produce the blue colour (BRAK 14) and ferrous (iron) to produce the translucent brown colour (BRAK 17).
The samples labelled BEY159 in Table 1 are unpublished glass mosaic tesserae found during excavations in Beirut, Lebanon by Dr Hans Curvers and dated to the first half of the 1 st century CE. All glasses are characterised by major levels of soda and silica; two contain significant levels of lead oxide (14.7 and 4.92). The tesserae contain much lower levels of magnesium and potassium oxides, indicating that the source of soda used was an evaporitic mineral instead of plant ash [21]. It is likely that the mineral used was natron which was used extensively to make glass by the Romans [22]. The colours of the opaque tesserae are caused by calcium antimonate crystals combined with copper to produce opaque turquoise (BEY159-T), lead antimonate crystals to produce opaque yellow (BEY159-Y) and lead antimonate combined with copper to produce opaque green (BEY159-G) [23].
A later Roman sample (denoted RAK2) is from the base of a 4 th -century aqua green goblet from the site of Rakkit in Palestine excavated by Professor Shimon Dar. The glass colour is typical for the period: thousands of blown aqua green blown glass vessels were made at the time in Palestine with Jalame being a key production site [24]. Like the tesserae from Beirut, this is a typical Roman natron glass with low levels of magnesium and potassium oxides associated with the use of the mineral natron used as a soda-rich flux. The green colour is due to a mixture of ferrous and ferric ions in the glass.

Computerised X-ray Tomography (CXT)
The samples were imaged using a High Resolution X-ray 3D Computed Tomography Microscope Instrument of model VeraXRM-510 (manufactured by Xradia Inc, Pleasanton, CA, USA). The field of view was 1024×1024 voxels per slice, with voxel resolution of 8 m.

Conventional gas sorption and overcondensation method development
The nitrogen adsorption/desorption isotherms were measured at liquid nitrogen temperature (77 K) using a Micromeritics ASAP 2020 physisorption analyser. Adsorption was measured over the relative pressure (P/Po) range of 0.01 to 0.995 for adsorption and desorption from 0.995 to 0.10 P/Po. The equilibration times used were 20 s and 60 s. Two different times were tested to ensure that the sample isotherm was fully equilibrated. The samples were outgassed at 140°C for 2 hours under vacuum prior to analysis.
The nitrogen overcondensation experiments were performed using a Micromeritics ASAP 2020 physisorption analyser using a method similar to that described by Murray et al. [25].
At the start of the overcondensation experiment, the pressure in the sample tube is increased to higher than the saturated vapour pressure of nitrogen. This pressure rise should facilitate sufficient condensation such that even the very largest pores are completely filled with liquid nitrogen at the start of the overcondensation desorption isotherm, which will almost certainly also involve some bulk condensation in the sample tube (see Figure 1). This bulk condensation is what is avoided in the conventional experiment. The required period to reach this stage is dependent upon the sample size and the pore volume. While it does not matter if the volume of condensate is much higher than that needed for complete pore filling (see Figure 1), the total duration of the experiment would be much longer in that case. Once complete pore-filling had been achieved, the pressure was lowered to just below the saturated vapour pressure of nitrogen such the bulk condensate vaporized completely while keeping all the sample internal porosity liquid-filled. Once this stage has been accomplished, the first data point on the overcondensation desorption isotherm can be measured. This point corresponds to the total pore volume of the sample. The pressure is then progressively lowered in small steps, and the rest of the desorption isotherm was obtained in the usual way.
Since the glasses studied here are low surface area materials, with relatively low nitrogen uptake, and only available in limited amounts (as is common for rare archaeological materials), it was necessary to develop the overcondensation method beyond previous work [8,25] to cope with such samples, and, thus, several variants of the overcondensation method were trialled to determine which was best for such materials. A key issue was to reduce the extraneous volume in the apparatus. A reduction in the free space, with a conventional sample tube, was achieved using a filler rod, or a filler rod plus 1-2 mm diameter ballotini, Unfortunately, these resulted in non-closure of the isotherm as shown in Figure 1. Using larger (4 mm diameter) ballotini helped the non-closure to disappear. However, the best approach found, and used subsequently (to obtain all overcondensation data in Results section), was to use a narrow cold trap tube just with a filler rod, which halved, or more, the warm and cold free space, compared to a conventional sample tube. The specific surface area was determined using the standard BET model which is only applied to the ISO recommended, more limited, fitting range of relative pressures from 0.05-0.3 [2]. refocusing pulse were measured on the glass specimens. The var delay was varied from 1×10 -experimentation, this is expected due to the disordered nature of local sodium chemical environments in glass and is supported by the observed dispersion of chemical shifts in onedimensional spin-echo experiments. The T2 determinations were obtained using fits of either a mono-exponential or bi-exponential decay function and the goodness of the fit was determined by regression analysis (R 2 , >0.99 in all cases), to establish whether a mono-or biexponential approximation was more appropriate. A normal distribution of the residuals was achieved during the analysis.

NMR
The 23 Na MAS and 3QMAS spectra (ν0 = 105.85 MHz) were achieved using a 100 kHz (2.5 µs, π/2) excitation pulse on Bruker 4 mm HXY probe achieving a MAS frequency of (νr =) 12-13 kHz, for all experiments. A solution reference of NaCl (aq, δiso = 0 ppm, 1 M) was employed to calibrate the pulse lengths and shift range. The T2 relaxation data was achieved using a 30 mm 23 Na MRI resonator, using a spin-echo sequence with a 10 kHz (25 µs, π/2) excitation pulse and a 5 kHz (50 µs, π) refocusing pulse, with a variable delay between the pulses. All data was processed using the Bruker Dynamics Centre and Topspin (Billerica, MA) software, and simulated using the DMFit [12].

SAXS
SAXS patterns of the powdered glass sample were measured on an Anton Paar SAXSPoint with a Cu K X-ray point source. A thin layer of the sample was placed on a piece of sticky tape and measured in a holder at 25 C under vacuum. A pattern from an empty piece of sticky tape was subtracted from the sample data to remove background scattering.  Table 1 that give the glasses their different colours. However, the morphology of these heterogeneities characteristically differs between the BEY159 and BRAK series. It is noted that, for BEY159-Y, the CXT slices reveal thin bands of both bright white and darker intensity within the grey matrix. The variations in X-ray absorbance in the bands are due to differences in the concentration of heavy elements present in the flux or colourant, for example [27]. BEY159-G and BEY159-T also show some evidence of fainter, solely darker or whiter banding, respectively, but these are much less pervasive and thinner than for BEY159-Y. This is because BEY159-Y has a lower melting temperature and a longer working period because it contains 14.7% PbO. This made the glass fluid for longer, and, in this case, strings of small lead antimonate crystals visible in the glass show how it has been stretched. In the CXT slices for BEY159-Y, it can also be seen that the black bubbles are often not circular, but frequently ovoid or lenticular in shape, with the long axis of the bubbles running parallel to the orientation of the bands of white and dark glass. One or two of the bubbles in BEY159-G also show some much smaller distortions from regular, circular shape. However, the bubbles in BEY159-G and BEY159-T are generally much more circular than those of BEY159-Y because the former contain PbO oxide levels of 4 and 0.02% respectively and the glass therefore has higher melting temperatures with shorter working periods. The banding and elongation of the bubbles in BEY159-Y is indicative of the direction of internal folding and stretching of the still molten, or softened, glass before it is cooled, and the patterns are then frozen in [26]. In contrast, the whiter and darker areas in the BRAK-series glasses are typically of the form of large, amorphous patches, or obviously follow cracks; in particular, BRAK 12 (Fig 2(d)) and BRAK 17 (Figs. 2(f) and A1(f)) reflect the presence of (darker) leached areas. In BRAK 14 (Figs 2(e) and A1(e)) the leaching follows cracks through the glass. This indicates that BRAK glasses were not poured in the way of the BEY159 glasses. The pouring lines observable in the BEY glass tesserae are due to their production process. Glass is first poured into a rectangular mould to the thickness of the glass tesserae and in so doing often produces lines of elongation in the tesserae produced.

CXT
Strips of the desired width of the tesserae are then cut from the block; individual tesserae are snapped off over a sharp edge.
The extensive network of cracks throughout the whole of the bodies of the chips of the BRAK14 and BRAK17 glasses may be indicative of brittleness resulting from insufficient annealing during fabrication [27]. A combination of large heterogeneities and numerous bubbles is often indicative of low-quality workmanship and low-quality glass [27]. In Figure   2(a), the CXT image of BEY159-G, one of the bubbles on the left-hand limb of the chip shows evidence of containing a particle within it. This is most likely to be a particle of polishing powder trapped in a sectioned bubble in the sample surface. The glass adjoining some of the bubbles and the cracks in BRAK14 and BRAK17 show less X-ray attenuation than the surrounding matrix. This is evidence for the leaching of elements from the glass surface, especially sodium, probably by water ingress in the burial environment along the cracks in the glass. This lack of durability is partly due to the relatively low levels of the network stabilizer, CaO, in the BRAK glasses (3.9%-5.8%) and also a range of other factors such as variable levels of humidity and temperature in the semi-desert burial environment of northern Syria.  Figure 3 shows the multilayer region of the experimental nitrogen adsorption isotherms obtained for chip samples of the BEY159-G. The equivalent plots for the rest of the BEYseries and the BRAK series of glasses are given in Figure A2 in the Appendix. components [4]. Figure 4 shows the multilayer region of the experimental nitrogen adsorption isotherm obtained for a fragmented sample of the BEY159-G with powder particle size less than 300 m and fits to the same type of models as used for chip data.   Tables 2, 3 and 4 show the fitted parameters for the standard BET, two-component fractal BET, and combined Henry's Law and Langmuir models for the glass adsorption data. Figure   5 shows a comparison of the respective specific surface areas obtained for each sample from fits of the three different isotherm models to the experimental data over relative pressure range 0-0.4.   From Figures 3 and A2, it can be seen that the standard BET model, fitted to the recommended range, always provides a poor fit beyond that range at high relative pressures, since it rapidly over-predicts the adsorbed amount. As mentioned above, the standard BET model neglects the impact of surface roughness and chemical heterogeneity. In contrast, the two-component fractal BET homotattic patch model, which accounts for both these aspects, fits the adsorption data for all the glasses well beyond the fitted relative pressure range. The two-component Henry's law and Langmuir homottatic patch model also follows the data beyond the fitted range, but, generally, not as far as the fractal BET model. As can be seen from Figure 5, in all cases, the two-component, homotattic patch models obtained larger specific surface areas than the standard BET method, and in some cases the area was twice that of the standard method. As can be seen in Tables 2-4, fragmentation of the BEY-159-G sample led to an increase in accessible surface area but did not significantly affect the parameters for the interaction strengths and fractions of each of the two-components of the homotattic models.

Gas sorption
In some cases the two-component fractal BET model fitted to relative pressures below 0.4 showed a positive deviation even at the highest relative pressures in the hysteresis loop region (see Figure A2(b),(d),(e)), which is impossible since that would imply multi-layer adsorption was greater than capillary condensation. Hence, in those case, the two-component fractal BET model was fitted to a larger relative pressure range up to a value of 0.9. Comparing the fractal BET model parameters between these two fits shows that the change in fitted range tended to increase the monolayer capacity (and thus surface area) and, generally, the fractal dimensions, but left the fractions of each component more or less the same. It is noted that, for the two-component fractal BET fits, for all samples C1<<C2, but for the BEY159 samples, D2 is generally greater than or equal to D1, whereas for the BRAK samples, D1 is generally greater than or equal to D2. From Equation (6), a larger value of apparent fractal dimension may be associated with a smaller size of individual patches (making up the overall fraction) of a given surface component.
The CXT data has shown that all the glass samples contain macroporous bubble-like pores, that appear isolated on the macroscopic scale (>10 m). Normally most macropores are beyond the detectability limit for conventional nitrogen gas sorption experiments, but they are potentially detectable using overcondensation experiments. Hence, those glasses showing a higher prevalence of macropores in CXT images were also studied using overcondensation.
The schematic diagram in Figure 6 shows the gas sorption isotherms anticipated from a conventional adsorption experiment for a system where the macroporous bubbles are connected to the exterior via nanopores. In the case of a conventional experiment, the isotherms would only detect sorption in the nanopores, as the macropores do not fill with condensate at all during the experiment. If an overcondensation experiment is performed, bulk condensation will be achieved (as indicated by the vertical solid line at the bulk condensation pressure in Figure 6), but if the macropores are disconnected they will still remain unfilled and undetected in the sorption isotherms (as shown in Figure 6(b)). If the macropores are connected to the exterior via nanopores, and an overcondensation experiment is conducted, then the macropores will fill at the bulk condensation pressure, and the step size in the boundary desorption isotherm will be much increased relative to the conventional experiment (as shown in Figure 6). The scenario shown in Figure 6(c) is what has been observed experimentally for shale rock samples with macropores shielded entirely by mesopores [9]. In the glasses, the scenario in Figure 6(c) would be expected if the nanoporosity had been formed by leaching of species from the glass matrix surrounding the bubble pores and cracks (as possibly suggested by CXT data given above) leading to a 'corona' of nanoporosity around the macroporosity.   Figure A3 in the Appendix.
The presence of hysteresis loops suggests the occurrence of capillary condensation, and thus the presence of some mesoporosity. For both BEY159-and BRAK-series samples, the hysteresis loop region is generally quite thin, except for BRAK14. Further, where conducted, the overcondensation full boundary curves are very similar to the desorption branches of the conventional isotherms, as would be expected for the scenario given in Figure 6(b). This suggests that there are no macropores that did not fill in the conventional adsorption experiment but which are shielded by smaller pores that condense at lower pressure than very close to the bulk step. Hence, the bubble pores evident in the CXT images are not connected to the exterior by the mesoporosity, otherwise a large increase in initial condensation, and a pronounced pore-blocking or cavitation step in the overcondensation boundary desorption isotherm, would be expected (as seen in Figure 6(c)). The small number of supermacroporous bubble pores with direct access to the exterior seen in the CXT images must fill and empty at pressures so close to the bulk that they contribute to the bulk step in the overcondensation experiment. Also shown in Figure 7 (and Figure A3) are the predictions for the position of the desorption boundary isotherm made using the data for the conventional adsorption isotherm and Equation (10), assuming exponents of 1.8 and 2.0 (corresponding to NLDFT and Cohan equations, respectively [16]), for single pore hysteresis. It was found that the single pore hysteresis equation makes good predictions for the data, with the exception of that for BRAK14, where the experimentally-observed hysteresis is much wider. This finding suggests the mesopore network in BRAK14 has substantial pore-blocking, suggesting significant (greater than a factor of 2 variation) corrugations in pore width. There may be pore width corrugations for the other samples but these must below the detectability limit (of a factor of 2 between neck and body sizes) imposed by the advanced condensation effect [16]. For the samples (besides BRAK14) where the hysteresis width was close to the single pore hysteresis predictions, the exponent of 1.8 gave a better fit to the BEY159 samples, whereas the exponent of 2 gave a better fit to the BRAK samples. This suggests the BRAK samples have a surface that is generally more wetting to nitrogen than the BEY159 samples.
In order to further study the connectivity of the macroporous bubble pores, for the sample BEY159-G, which had the highest prevalence of these pores in the CXT images, conventional and overcondensation nitrogen sorption isotherms were obtained following fragmentation to a powder particle size <300 m, which is less than the typical distance between many bubble pores in the CXT images, and these are shown in Figure 8. From Figure 8, it can be seen that the overcondensation boundary desorption isotherm is still very similar to the conventional desorption isotherm, but the total uptake at the top of the boundary isotherm has increased slightly from ~8-9 to ~12 cm 3 (STP)g -1 . This finding suggests that whereas decreasing particle size has made some mesopores accessible, this is not the case for any further large macropores. Again, any bubble macropores with access to the exterior must fill and empty at, or very close to, bulk condensation. Since the majority of the conventional and overcondensation desorption isotherms suggest that the hysteresis in the gas adsorption data has a single-pore origin, then it is appropriate to use the Cohan equation for desorption to obtain the pore size distribution. The Barrett-Joyner-Halenda (BJH) desorption pore size distributions are shown in Figure 9.

NMR
Sodium-23 NMR spin-echo T2 relaxometry experiments were performed on the BEY159, BRAK and RAK series of glasses, and the results are shown in Figure 11, along with fits of mono-exponential and bi-exponential models for the signal decay. The relaxation time parameters obtained from these fits are shown in Table 5. From Figure 11 and Table 5, it can be seen that the relaxation time data for the BEY159 series samples give the best fit to a biexponential model, whereas the BRAK series and RAK2 samples best fit to monoexponential models.  Table 5. Fitted parameters for single and double exponential decay models for spin-spin (T2) relaxation time data for sodium-23 in glass samples shown in Figure 11.   Figure 12 shows the 3QMAS data only for selected glass samples (as the technique is destructive, in contrast to non-destructive relaxometry). The data for the RAK2 and BRAK samples suggest a single resonance, while that for BEY159-G shows two resonances, centred at MAS of -2 and -9 ppm, with the latter broader than the former. However, the two resonances overlap (likely chemical environments as sodium is a counter-ion in glasses). The presence of the two components, for BEY159-G is consistent with two sodium environments giving rise to two different relaxation times in the T2 relaxometry data in Figure 11. The BEY159 sample has quite narrow NMR lines, suggesting a more ordered material. The lines for the other two types of glass are broader, as expected with a glass. Figure 13 shows typical FEG SEM images for the BEY159-G, BEY159-Y, and BRAK12

Electron microscopy
samples. Further examples of FEG SEM images of these samples are given in Figure A4 in the Appendix. Nanopores of a size similar to that from the BJH desorption PSDs can be seen in the images as dark ovoids.

DISCUSSION
It has been seen, from the CXT and gas sorption data, that the BEY159-series and BRAKseries glass samples possess both macroscopic, bubble-like pores, and nanoporosity. The form of the sorption hysteresis loops is that expected if the hysteresis were single pore in origin, which is typically associated with open cylindrical pores. In contrast, the SEM data suggests that the nanopores form a more irregular network. However, the hysteresis can still appear to be single pore origin if the size variation within the network is only of the order of a factor of two (as it mostly is according to PSDs) since then shielded pore bodies would fill at the same time as the shielding neck via advanced adsorption, and empty at the same time due to pore-blocking, and, thus, the isotherm is essentially just that of the neck alone [16].
The overcondensation experiments suggested that the two levels in the hierarchy of porosity are not connected, and the macroporous bubbles are not like 'vugs' within a sea of nanopores (as often occurs in rocks [28]), and, thus, the latter do not form a 'corona' within the leached matrix around the bubbles. This is because, if the macropores visible in the CXT images were shielded but inter-connected by nanopores, then the overcondensation experiments would have enabled them to be detected in the overcondensation desorption isotherm, as a substantial growth in the capillary evaporation in the mesopore range relative to the conventional isotherms, as has been observed previously for large macropores shielded by mesoporosity in shales [9]. The fragmentation of the BEY159-G sample led to an increase in the accessible nanoporosity and, (since, the hysteresis remained of single pore origin) no detectable increase in macroporosity. This suggests the bubble pores exposed on the surface of the sample particles (as seen in Figures 2(a),(b),(c)) fill and empty with the bulk condensate.
It has been shown that the conventional BET analysis of gas adsorption data, according to the ISO method, does not fit the data outside the fitted range, suggesting the model is inadequate.
It was found that both two-component Henry's law and Langmuir isotherm, and twocomponent fractal BET, homotattic patch models, were better able to follow the multilayer adsorption up towards the region of capillary condensation even when only actually fitted to a much lower relative pressure range. This suggested that this two-component model was better than the standard BET approach, which implies that the glass surfaces possess significant chemical and/or geometrical heterogeneity sufficient to severely affect the accuracy of pore space descriptors. Several aspects of the gas sorption data show differences between the BEY159 series and BRAK series samples.  Figure 14 suggest there is a strong correlation between the overall iron content and the measured overall average T2. No such degree of correlation was found for the other elements present, even the paramagnetic species. The iron present in the glass samples will have originated from minerals within the sand used to make the glass, and will have, subsequently, permeated throughout the glass when it initially became molten leading to intimate mingling of iron and sodium-rich components, such that the (average) sodium T2 was determined by the overall iron content.
However, the two-component fit needed for the BEY159-series samples suggests that the two different chemical environments of the sodium are associated with slightly different local concentrations of iron.  Table 1) for each glass. The lines shown are straight-line fits to data shown. Figure 15 shows the variation with overall iron content (from Table 1) of the weighted average heat of adsorption (calculated using equation (3) The above findings from the homotattic patch model can be explained as follows. The nanopore surface for all types of glass studied is considered to have a distribution of high energy adsorption sites (associated with the Langmuir patches), and some low energy sites (associated with some of the Henry's law patches), as shown schematically in Figure 16. In each type of glass, the available iron is considered to occupy, or block access to, the very highest energy sites first, and, then, progressively block lesser high energy sites until the iron available at the pore surface of a particular type of glass is exhausted (which is proportional to overall iron content). For each type of glass, this disposition of iron across the surface will leave a rump of remaining higher energy (Langmuir fraction) sites, with the numbers left available depending upon the iron content of the particular type of glass. Nitrogen gas will interact more strongly with the remaining higher energy sites (which form the Langmuir patch) and the sites blocked by iron will be perceived by nitrogen as part of the Henry's law component. This model explains how the average heat of adsorption of nitrogen and B are inversely correlated with the iron content. Figure 15. Variation with overall iron content (from Table 1) of the weighted average heat of adsorption (calculated using equation (3) Figure 16 was found to be 0.9987, which is significant at the 95% confidence level, even for such a small sample. This suggests that the fractions of these two phases in each model for BEY159 glass are related, thereby suggesting the respective model parameters have physical meaning. Given the more general correlations found with iron content above, the phases may relate to the relative local distribution of iron between bulk matrix (probed by NMR) and pore surface (probed by gas sorption) for the BEY159-series glasses. Gas overcondensation allows a check that larger mesopores and macropores are not missed from the conventional gas sorption pore size distribution. It is thus noted, from a comparison of the ultimate mesopore volumes of the glasses shown in Figure 9 with the elemental compositions in Table 1, that the glasses with the highest alkalinity (sodium plus potassium), BRAK-12 (20.05 %) and BRAK-14 (20.58 %), have the highest mesopore volumes, while the glass sample with the lowest alkalinity, BEY159-Y (12.91 %) has the lowest mesopore volume. This is consistent with the previous suggestion [29] that the variation in the nanoscale porosity in ancient glasses relates to alkali content. It is noted, from the CXT images of BRAK12 and BEY159-Y in Figure 2, that the glass chips with the biggest difference in mesoporosity are of similar (macroscopic) size, and neither shows evidence of highly pervasive macroporous cracks. Hence, their difference in mesoporosity is unlikely to be due to increased accessibility due to simply particle size differences, or weathering. The EPMA analyses were for a surface analysis of a small unweathered part of the glass. This also suggests the heterogeneities in iron content found above do not result from leaching due to weathering.

CONCLUSIONS
The standard BET method has been shown to be completely unsuitable for providing accurate surface areas for making comparisons between samples, since it does not account for the degree of surface heterogeneity seen for ancient glasses. In contrast, the two-component homotattic patch model not only provided a more satisfactory fit to the gas sorption data, but the model parameters have also been demonstrated to possess independent physical meaning, since the weighted-average heat of adsorption correlated with the overall iron content of the glass. The iron content of the glasses was also shown to predominantly determine the average sodium-23 relaxation time. The type of glass (BEY-series) that had the most clear suggestion of two different sodium environments from the NMR data, also had the best correlation between the fractions of the respective NMR and gas sorption phases. The Bronze Age and BEY-series ancient glass samples have also been shown, by CXT and FEG-SEM, to possess a hierarchy in porosity over length-scales of both ~100s microns, and a few nanometres. A new, non-destructive, gas overcondensation characterisation method suitable for low surface area materials was developed to study the sample-wide inter-relationship, between the various levels of the multi-scale porosity, that microscopy and imaging cannot provide. This new method was used to show that the large bubble macropores were isolated from the nanoporous network in all cases, and, thus, the latter did not represent the hypothesised 'corona' around the bubbles and cracks formed by leaching, but it had another source possibly related to sample alkalinity.