Diffusion of Hydroxylated Si Vacancies in Olivine, and Its Relevance to Si Diffusion

The diffusion coefficient associated with the 3,612 cm−1 infrared absorbance band in H‐bearing olivine, attributed to a defect involving a vacant Si site balanced by four H+ (i.e., (4H)Si× ${(4\mathrm{H})}_{\mathrm{S}\mathrm{i}}^{\times }$ ), was determined in H in‐diffusion experiments at 1200°C–1400°C and 1.5–3.0 GPa in piston‐cylinder and multi‐anvil apparatuses. The retrieved values are in good agreement with experimentally determined diffusion coefficients ascribed to the same process from previous experiments, despite those being conducted at atmospheric pressure and representing out‐diffusion rather than in‐diffusion. Overall, the diffusivity of the Si vacancy is several orders of magnitude lower than that of the M‐site vacancy (∼4 orders of magnitude at 1000°C), with a higher activation energy (∼450–500 vs. 200–250 kJ mol−1), when both are charge‐compensated solely by H+. From these data, the diffusion coefficient of Si in H‐bearing olivine can be estimated using the relationship between vacancy concentration, vacancy diffusivity, and diffusivity of an ion moving by a vacancy mechanism. However, these diffusivities are inconsistent with published results for Si diffusivity, with a considerable discrepancy of around two orders of magnitude. This suggests that simple relationships between diffusion coefficients of vacancies and species moving by vacancy diffusion mechanisms may be inappropriate, highlighting again the complexity of H incorporation and diffusion mechanisms in olivine.


Introduction
Largely owing to their high charge, silicon vacancies in forsterite, in the pure Mg-Si-O system, are energetically unfavorable (Brodholt & Refson, 2000;Muir et al., 2022;Smyth & Stocker, 1975;Stocker & Smyth, 1978). This suggests that it is unlikely that Si diffusion in pure forsterite (Mg 2 SiO 4 ) occurs by a simple vacancy mechanism; rather, a mechanism involving interstitial sites has been proposed (Bejina et al., 1999;Brodholt & Refson, 2000;Houlier et al., 1990). Understanding the rates and mechanisms of Si diffusion is important for considerations of deformation, as Si diffusion may be a controlling factor for creep at certain conditions (Boioli et al., 2015;Gasc et al., 2019;Houlier et al., 1990;Jaoul et al., 1981), although there is a long-standing discussion regarding specific mechanisms, and the relative importance of Si versus O (Hirth & Kohlstedt, 2004;Kohlstedt & Ricoult, 1984;Ryerson et al., 1989). However, when H is added to the pure Mg-Si-O system, the Si vacancy can be stabilized by full hydroxylation, that is, the formation of four OH − groups by adding H + to the O 2− surrounding the vacant tetrahedral site (Balan et al., 2011;Braithwaite et al., 2003;Lemaire et al., 2004;Muir et al., 2022;Qin et al., 2018;Walker et al., 2007). This forms charge-neutral defect Mg 2 H 4 O 4 /(4H) × Si in endmember and Kröger-Vink notation, respectively (a.k.a. the "hydrogarnet" defect). A more complete description of the defect in Kröger-Vink notation is Si is used herein for brevity. In Fourier transform infrared (FTIR) spectroscopy, this defect is expressed as a series of bands in the ∼3,500-3,630 cm −1 wavenumber range, with the strongest generally being a band at 3,612 cm −1 (Figure 1, Berry et al., 2005;Lemaire et al., 2004;Matveev et al., 2001).

Where
(4H) × Si defects exist, exchange between these defects and Si 4+ should enable a vacancy mechanism to occur for Si diffusion, which, according to calculated energetics of the H-free Si vacancy, can be effectively considered as impossible in the H-free system. If this vacancy mechanism is possible, the diffusivity of Si will be a function of both the concentration and diffusivity of the (4H) × Si defect. Increasing either or both of these will increase the probability of Si × Si ↔ (4H) × Si exchange. This relationship, is consistent (in terms of direction but not magnitude) with the data set of Fei et al. (2013), and the general increase in Si diffusivity between H-free (Bejina et al., 1999;Dohmen et al., 2002;Fei et al., 2012) and H-bearing (Costa & Chakraborty, 2008;Fei et al., 2013) Si diffusion experiments, discussed in more detail below. Whilst there are several studies describing the effect of temperature (T), pressure (P), silica activity ( SiO 2 ) water fugacity ( H 2 O ) and oxygen fugacity ( O 2 ) on the concentration of the (4H) × Si defect, or at least the absorbance of the 3,612 cm −1 band even if interpreted differently 2 of 16 (Bai & Kohlstedt, 1993;Berry et al., 2005;Lemaire et al., 2004;Matveev et al., 2001;Mosenfelder, Deligne, et al., 2006;Tollan et al., 2017;Withers & Hirschmann, 2008), only one study describes the diffusion of this defect (Padrón-Navarta et al., 2014). Given the potential importance of this defect for Si diffusion and hence deformation rates and mechanisms, it is imperative that these published diffusion data are independently validated.
In this study, the problem of (4H) × Si diffusion is revisited, with H in-diffusion experiments conducted between 1200°C and 1400°C at 1-3 GPa. Despite considerable differences in experimental (in-diffusion at high pressure vs. out-diffusion at atmospheric pressure) and analytical protocols (resolving (4H) × Si vs. distance profiles rather than bulk (4H) × Si loss as a function of time), the data from this study agree well with the Padrón-Navarta et al. (2014) study, to the extent that the data can be fitted together to give a generally applicable Arrhenius relationship. Attempts to link the measured diffusion coefficients with those of Si are unsuccessful, however, which suggests that simple models relating the diffusivity of the Si vacancies with Si diffusion may be inappropriate, given the complexity of H transport and incorporation mechanisms in olivine.

Experimental
In a perfect experiment, the aims would be to (a) partially hydroxylate single crystals of synthetic forsterite or natural olivine at a given and known P, T, SiO 2 , H 2 O , and O 2 ; (b) form a H diffusion profile that is sufficiently long to measure by FTIR spectroscopy but not so long that concentration plateaus are removed in crystal cores; and (c) recover the single crystal in an acceptable state for FTIR spectroscopy. It was not generally possible to accomplish all of these aims simultaneously at the P-T-t conditions required for this study, such that , are unknown at run conditions. This then means that it is not possible to determine O 2 using the presence of graphite (used as a capsule liner) and H 2 O. Based on these initial findings, the following method is the result of several iterations and compromises. All experiments and analyses were conducted at the Lamont-Doherty Earth Observatory (Columbia University) in New York, USA.
Experimental conditions are detailed in Table 1. Single crystals of synthetic forsterite (Czochralski-grown) or San Carlos olivine were oriented, then cut into cubes of ∼1.2 mm on each side, with the faces of the cubes corresponding to the (001), (010), or (100) planes. These were oriented using a combination of planes present on the as-grown material and Si-O overtones in polarized FTIR spectra. The forsterite was from the same stock used by Zhukova et al. (2014)-trace element compositions are presented in that study. The starting forsterite can be considered as absolutely H-free, given its growth method, and indeed no absorbance bands were visible in the OH-stretching region in FTIR spectra. Similarly, the San Carlos olivine starting material showed no absorbance associated with OH, consistent with the starting materials used by Tollan et al. (2018).
a Experiment done using a multi-anvil apparatus, all others were used piston cylinder. Column "Powder": per = periclase, brc = brucite, fo = forsterite, and tlc = talc. b A band present as a shoulder on another, so is only identified by resolving spectra into Gaussian peaks. c The band that is likely related to spectral contamination from brucite present at the crystal edge.

of 16
(1200°C-1400°C, monitored using a type "D" W-Re thermocouple), with the pressure gradually increased to the final pressure (1-2 GPa) during the temperature up-ramp. Following an isothermal, isobaric dwell, the temperature was ramped down at 100°C/min to room temperature, then the pressure was released over at least 10 min. One experiment was also performed at 3.0 GPa in a Walker-type multi-anvil apparatus in an MgO assembly with a fecralloy (Fe-Cr-Al alloy) heater. A miniaturized version of the same capsule design was used as for the piston-cylinder experiments.
The whole capsules were then mounted in epoxy resin in 2.5 cm discs, and ground down from one end until the (010) face of the crystal was just exposed. Then, another ∼400 μm was removed parallel to (010), and the face of the epoxy mount polished down to 1 μm Al 2 O 3 , then the reverse side of mount was ground down to around 400-700 μm thickness, and polished again, thus yielding a doubly polished thick section ( Figure 2b) representing approximately the core of the original 1.2 mm cube.

Analytical
Fourier transform infrared spectra were recorded using a Thermo Nicolet i10 microscope. All spectra used for the determination of diffusion coefficients were recorded using unpolarized light, with polarized spectra only recorded for the purpose of verifying the orientation of the crystals post-experiment. This was done by comparing the Si-O overtone region of polarized spectra with the reference spectra of . Spectra were baseline corrected using a routine implemented in MATLAB, where a parabolic (U-shaped) function was added to each spectrum, then a bounding curve was drawn around the base of the modified spectrum, smoothed, then subtracted from the spectrum. This method removes broad oscillations, which were commonly observed in the spectra (likely due to internal fractures) but preserves sharp peaks. The method was tuned such that the relatively broad bands at around 3,150-3,220 cm −1 , present in some spectra and assigned to a hydroxylated M-site vacancy (e.g., Berry et al., 2005;Lemaire et al., 2004), were not inadvertently deleted by baseline correction. In fact, in this study, the choice of baseline correction routine has a negligible effect on the results, as the band of interest is the relatively tall, sharp band at 3,612 cm −1 , and it is not necessary to convert integrated absorbance to H 2 O content for the purpose of the simple diffusion modeling employed herein. Spectra were then corrected to 1 cm thickness using thickness measured with a digital micrometer, and multiplied by 2 to give an approximation of total absorbance. The "2" factor is an approximation (likely a slight underestimate) and is only valid for unpolarized spectra recorded from a thick section of an olivine crystal prepared parallel to its (010) plane, given that the majority of absorbance associated with the O-H bonds of interest would be observed using polarized light when the electric vector (E) is parallel to [001] or [100] (e.g., Jollands et al., 2021;Lemaire et al., 2004;Mosenfelder, Sharp, et al., 2006;Withers et al., 2012).
In order to extract the absorbance of the 3,612 cm −1 band, the 3,400-3,750 cm −1 region of the spectra was resolved into a series of Gaussian peaks. This is generally necessary given that the band often overlaps with adjacent bands, so simple numerical integration may be inappropriate.

Approximating the Effective Beam Size
FTIR spectra were recorded using a 30 × 150 μm aperture, with the long axis parallel to the crystal edge at the start/end of the profile. Rectangular apertures allow simultaneously high spatial resolution, from the short axis, and generally increased quality of spectra (relative to the equivalent square) from the long axis. The actual size of the beam, in terms of the area represented by a single spectrum, is likely to be somewhat larger than this. This is largely due to the divergence of the beam within the single crystal, with the magnitude of this effect increasing with increasing sample thickness. In order to be able to determine the extent of this convolution, a single crystal of olivine from Pakistan (as described by Gose et al. (2010)) was cut into slabs perpendicular to [001], then each slab was cut in half. Each piece was then polished on one (001) face. The polished (001) faces from the two halves of each slab were pressed together with one rotated 90° relative to the other around the [001] axis, and epoxy was painted around the sides of the couple, thus joining the two pieces together. Each of the couples was then cut into several slices, which were then mounted in epoxy resin in 2.5 cm discs, with the (001) faces perpendicular to the face of the disc. These were then prepared as doubly polished sections of various thicknesses: 223-245, 321-330, 386-396, 477-486, 547-559, 638-647, and 849-866 μm, with the ranges representing at least five measurements across the sections. The aim of this preparation is to obtain a couple with a perfect step function in absorbance in FTIR spectra, which can be profiled across to determine the convolution effect associated with the diffusion profile analyses. Whilst the crystal on either side of the interface is nominally identical, the 90° rotation coupled with the anisotropy of OH bonds in olivine means that the FTIR spectra on either side will be distinct.
A series of unpolarized spectra were then recorded in a transect across the couple, as above, with the long axis of the aperture (30 × 150 μm) oriented parallel to the interface. Examples are shown in Figure 3. The band centered at 3,612 cm −1 showed the greatest difference in absorbance across the couple, so the absorbance at this wavenumber was extracted. The absorbance-distance curve was then fitted to the cumulative distribution function of the normal distribution, assuming that the beam-sample interaction is approximately Gaussian, to extract the standard deviation, and then the full width at half maximum (FWHM), associated with the Gaussian curve.
A second test was conducted to directly simulate the experiments. One piece of pure forsterite was placed into a pellet press along with brucite powder and a few drops of polyvinyl alcohol glue. This was then pressed and dried, then mounted in epoxy, and prepared into a doubly polished section (633 μm thick), as above, thus the section contained a single crystal of forsterite surrounded by pressed brucite. A series of spectra were then recorded in a linear transect across the crystal, as above. Unlike the tests described above where the spectra on either side of the couple are relatively similar, the aim of this test was to determine the distance over which the signal from adjacent brucite could be observed in dry brucite-free forsterite.
The FWHM of the interaction volume is consistently greater than the nominal aperture width (30 μm), although at ∼230 μm sample thickness, the FWHM approaches the 30 μm limit. The relationship between sample thickness and FWHM is broadly linear over the studied thickness range. Similar behavior is observed for the olivine-olivine couple and the pressed brucite powder-forsterite couple, albeit with only one profile of spectra recorded from the latter.
The reason for the scatter in Figure 4 is not clear, but may relate to linear transects being recorded at different positions of the crystal, or the focal plane being in different positions in the vertical axis. Because the method involved recording FTIR spectra, thinning the samples, recording spectra along new profiles, thinning, and so on, it is not possible to repeat individual measurements.
For the purpose of assessing convolution effects in this study, the maximum and minimum FWHM were simply estimated by reading off the graph in Figure 4 (e.g., 50-75 μm FWHM at 500 μm thickness), with these values used in determining the minimum and maximum diffusion coefficients, respectively. A future systematic study of this effect would be useful, especially including considerations of other analytical parameters (e.g., numerical aperture, physical aperture, focal position, material being analyzed, hardware). The FWHM values obtained here are similar to the 30 μm FWHM determined by Ni and Zhang (2008)   . Examples of spectra recorded along a linear transect across an olivine-olivine couple, ∼550 μm thickness, with a cartoon in the inset. The olivine crystals on either side are nominally identical, but one is rotated 90° relative to the other, hence the difference in absorbance between bands on either side of the boundary. Spectra (i) and (iii) represent the end members, and (ii) is a mixed spectrum recorded near the boundary. The difference in baseline shape likely relates to some contamination from organics/polishing compound in the small crack between crystals. 6 of 16

Extracting Diffusion Coefficients
Profiles of integrated absorbance of the 3,612 cm −1 band as a function of distance were fitted to one of two equations to derive band-specific diffusion coefficients, using unweighted nonlinear least squares regression. It was necessary to apply two diffusion coefficients (Ds) given profile shapes. This is discussed in detail below but, basically, profiles comprised a steep section (slow diffusion) and a shallow or flat section (fast diffusion). For full (rim-to-rim) profiles, Equation 1 was used: or Equation 2 for half (core-to-rim) profiles: Where: t is time (s); 2l is the length of the crystal (m); x is distance (with 0 < x < 2l in Equation 1 and 0 < x < ∼l in Equation 2); * slow is the diffusion coefficient associated with a slow mechanism; * fast is the diffusion coefficient associated with a fast mechanism. A is the integrated absorbance (cm −2 ), A rim,slow and A rim,fast are the absorbance associated with the slow and fast mechanisms projected to the crystal edge, respectively. Several different D notations are used in this manuscript: • * slow : D used to fit the slow mechanism, without any correction • * fast : D used to fit the fast mechanism, without any correction • D slow : * slow , corrected for beam convolution effects Equation 1 represents the superposition of modified versions of Equations 4.19 and 2.14 of Crank (1975), the latter adjusted to represent diffusion from two parallel faces of the same slab, as in, for example, Demouchy and 7 of 16 Mackwell (2006). Equation 2 represents two instances of Equation 2.14 of Crank (1975) with different diffusion coefficients superimposed. These equations are simplifications when considering the fast diffusivity * fast given that (a) it is one dimensional and the high diffusivities associated with * fast are such that a 3D model would be more appropriate and (b) the fast diffusion would be better described using a diffusion-reaction model-this is discussed below. It is for these reasons that no * fast values are presented- * fast is used only for the purpose of curve fitting. These simplifications are not problematic for determining * slow . In one experiment (FS8, 1400°C), there was a decrease in the absorbance of the 3,612 cm −1 band toward the crystal rims superimposed upon the increase from rim to core expected for in-diffusion, forming M-shaped profiles. Non-constant boundary conditions cannot be ruled out given the experimental design, so these areas of the profiles were not included in model fits-only the central U-shaped section of the curve was fitted (the profile and fit are presented in Figure S1).
Given that the beam-sample interaction produces error function forms when analyzing over a known step function, the convolution effect can be accounted for simply using the Ganguly et al. (1988) or Arnould and Hild (2002) formulations, or a numerical method such as that incorporated into Diffuser software (Wu et al., 2022). These were both determined for resolving beam-sample interaction when using electron microbeam techniques, but are entirely appropriate here. Therefore, the * slow values were corrected using the simple relationship: where σ is the FWHM defined above (Figure 4) divided by 2.355.

Uncertainties
Assuming that * slow values can be derived by considering 1D diffusion, the main sources of uncertainty associated with extracting slow (convolution effects removed) from a given profile are, in no specific order: 1. The uncertainty associated with the curve fit; 2. Uncertainty in the position of each measurement relative to the crystal edge as a result of crystal faces not being exactly perpendicular/parallel to one another/the face of the epoxy disc; 3. Profiles being artificially lengthened by the beam-sample interaction volume being non-zero.
To account for (a), the uncertainties (95% confidence intervals) associated with curve fitting were derived using MATLAB function nlparci, using the covariance matrix determined during non-linear least squares regression by nlinfit. For (b), the distance from the crystal edge for each spectrum was determined twice, one assuming that the first position was 20 μm from the edge, and once assuming 70 μm. These values are related to the edges of the crystal not being exactly perpendicular to the face of the thick section, and the 20-70 range represents the potential range of offsets for all profiles, determined using a petrographical microscope by focusing on both the top and bottom of the section. The data were fitted twice, using these two distance vectors, with which * slow and the 95% confidence intervals were determined. For (c), the maximum and minimum FWHM associated with beam-sample interaction was estimated using the results of the tests described above (section "Beam-sample interaction volumes"). The uncertainties associated with the three sources were then summed. The uppermost value of the 95% confidence interval for * slow determined using the highest offset was corrected for beam convolution using the minimum estimated FWHM. Likewise, the lowermost value of the 95% confidence interval for * slow determined using the lowest offset was corrected for beam convolution using the maximum estimated FWHM. For simplicity, the best fit log 1 0 slow is then assumed to be equidistant from these two points.

FTIR Spectra and Resolved Bands
The FTIR spectra recorded from crystals post-experiment show a plethora of bands (Table 1, Figure 5). It is assumed in the following that the only elements present at any non-negligble concentrations in the pure forsterite are Mg, Si, O, and Al (the latter at levels of tens of wt. ppm, Zhukova et al., 2014), and that their valence states are always 2+, 4+, 2− and 3+, respectively.

of 16
The only bands that are clearly visible in all recorded near-interface spectra are those at 3,612 and 3,624 cm −1 , with the latter consistently smaller than the former. All near-interface spectra from experiments using forsterite also show a triplet, with maxima at 3,567, 3,578, and 3,591 cm −1 , along with shoulders associated with bands centered at 3,571, 3,585, and 3,597. Whether this triplet is present in spectra from the experiment using San Carlos olivine is not clear, given the overlap with the high-absorbance bands at 3,534 and 3,572 cm −1 , attributed to the "Ti-clinohumite" type defect, Berry et al., 2005;Tollan et al., 2017). The bands at 3,567, 3,578, 3,597, and 3,612 cm −1 can be confidently assigned to (4H) × Si based on forsterite synthesis experiments conducted under different SiO 2 conditions (Lemaire et al., 2004;Padrón-Navarta et al., 2014). Specifically, these bands all show considerably stronger absorbance in forsterite crystals synthesized at low SiO 2 , which is consistent with a defect involving a Si-site vacancy. In the same wavenumber region, there are also resolvable bands at 3,571, 3,585, 3,591, and 3,624 cm −1 , which, based on their wavenumbers, could also be associated with the Several spectra also show broad bands at 3,160 and 3,218 cm −1 , which are generally attributed to fully hydroxylated M-site vacancies, (2H) × Mg as described above. Whilst the bands show high absorbance in some spectra, 9 of 16 such as the high SiO 2 experiment FS15, the contribution to the total OH content associated with this defect is probably minor. This is because the absorption coefficient for converting absorbance to concentration via the Beer-Lambert law is likely to be wavenumber/defect dependent, with higher absorption coefficients at lower wavenumbers (Balan et al., 2011;Kovacs et al., 2010;Libowitzky & Rossman, 1997;Paterson, 1982).
Bands at 3,321, 3,344, and 3,350 cm −1 can be attributed to a defect where Al 3+ replaces Mg 2+ , with the charge accommodated by a vacant M-site and a H + , that is, Berry et al., 2007). The small bands between 3,382 and 3,529 cm −1 have not been conclusively identified, but may relate to some combination of Al 3+ (or other trivalent cation), H + and Si sites, for example, (Jollands et al., 2021). A broad band centered at around 3,700 cm −1 is present in the near-rim spectra of many experimental charges. This is attributed to brucite, which is present in the buffer powder surrounding the crystals.
These band assignments are useful for qualitatively describing diffusion mechanisms, below, but for the purpose of determining diffusion coefficients related to the hydroxylated Si vacancy (4H) × Si , the 3,612 cm −1 band is sufficient. It is also the most appropriate for comparing different experiments, as it is the main band that is clearly present in all experiments. Additionally, as it is generally the strongest band in the ∼3,600 cm −1 region, it is least affected by overlap with other bands, and potential differences in baseline corrections will have only minor effects on the determined diffusion coefficients. In any case, fully resolving profiles, described below, shows that the bands assigned to (4H) × Si all appear to have similar effective diffusion lengths. Descriptions of diffusion coefficients herein generally refer to this band only.
One important note is that, according to the spectra in Figure 5, there is no apparent relationship between aSiO 2 , temperature or pressure on the concentration of (4H) × Si . This is very likely due to these experiments being unbuffered in terms of fH 2 O, as described in the methods.

Profile Shapes
Profiles of absorbance of the 3,612 cm −1 band versus distance from the crystal edge show two common features in all experiments ( Figure 6). The first is core concentrations that are non-zero, and the second is a relatively sharp decrease in 3,612 cm −1 absorbance within the first few hundred μm of the crystal edge, followed by either a gentle decrease in concentration toward the crystal core, or a flat profile ( Figure 6). This is the reason for fitting using equations which effectively superimpose two solutions to the diffusion equation. In some cases, it is notable that even this method does not fully describe the data. In Figure 6b, for example, the model appears to underpredict and overpredict at around 250 and 350 μm from the edge, respectively. This might indicate a process involving a combination of diffusion and reaction, which is expected to produce such deviations . . In this case, the concentration profile in the core region is effectively flat, therefore the * fast value is indistinguishable from ∞, that is, meaningless. This does not affect the * slow value.
Whilst it is not relevant to the primary aim of this study, it should be noted that resolving profiles into a series of Gaussian bands (Figure 7) shows that profiles associated with different bands have different forms, both in terms of the length scales over which concentrations increase/decrease, and the relative changes in absorbance. One example is shown in Figure 7 to demonstrate the complexity revealed by resolving spectra into individual bands.

Diffusion Coefficients
Diffusivity associated with the 3,612 cm −1 band ( Figure 8) increases with temperature, as expected, but there is no resolvable effect of pressure between 1.5 and 3.0 GPa. This is not a perfect comparison, however, given that the experiments were conducted in piston cylinder and multi-anvil apparatuses, respectively. The diffusivity is the same within uncertainty in San Carlos olivine and pure forsterite, and whether the SiO 2 is buffered by periclase (per, brc powder) or enstatite (tlc, fo source Taking only the D slow ||[001] values, fitting an isobaric Arrhenius relationship (assuming no pressure effect), log 10 D slow = log 10 D 0,slow − (E A /(2.303RT)) yields E A (activation energy) of 489 ± 37 kJ mol −1 , and log 10 D 0,slow (pre-exponential factor) of 4.6 ± 1.3 m 2 s −1 . The quoted uncertainties are approximations of 2σ, based on the assumption that the uncertainties for each fit, described above, also represent approximately 2σ. Because the diffusion coefficients determined in this study appear to align with those of Padrón-Navarta et al. (2014), a fit to all data from this study (both orientations) and the relevant data from Padrón-Navarta et al. (2014) yields an E A of 490 kJ mol −1 , and log 10 D 0,slow of 4.6 m 2 s −1 . The values hardly change from the D slow ||[001] fit, given that the relatively large uncertainties quoted by Padrón-Navarta et al. (2014) mean that these data have little effect on the overall weighted fit. In this case, formal uncertainties are not determined given the combination of different diffusion directions and the difference in experimental methods between studies, but a reasonable uncertainty estimate appears to be ±0.3 log units around the best fit between 1000°C and 1400°C.

Incorporation and Diffusion Mechanisms
All profiles of (4H) × Si versus distance showed forms indicative of fast and slow diffusion superimposed. This is evident either as relatively steep drops in concentration followed by shallower decreases (Figure 6a) or elevated background concentrations (Figure 6b). Elevated background concentrations resulting from fast diffusion in the third dimension, rather than a fast pathway, can be ruled out-this would require diffusion parallel to the [010] direction to be orders of magnitude faster than along [001] or [100], which has not been observed for any diffusion in olivine. In general, cation diffusion is fastest parallel to [001], [100] or close to anisotropic (e.g., Coogan et al., 2005;Demouchy & Mackwell, 2003;Dohmen et al., 2007;Kohlstedt & Mackwell, 1998;Petry et al., 2004;Zhukova et al., 2014).
The diffusion coefficients associated with the slow mechanism, * slow , are apparently insensitive to SiO 2 and olivine composition. This would be the  show in-diffusion but with different relative intensities. The unidentified band centered at 3,641 cm −1 shows apparent out-diffusion. These curves have been smoothed using a five-point moving mean; hence, there is a slight difference between these data and the data in Figure 6. 11 of 16 expected outcome for diffusion associated with exchange between (4H) × Si and Si × Si , that is, the movement of Si 4+ into a nearby vacant site, along with the movement of 4H + in the opposite direction. This is despite there being a known effect of SiO 2 on (4H) × Si concentration (e.g., Lemaire et al., 2004), which suggests that (4H) × Si diffusion is concentration-independent. Note that this should not be taken as a suggestion that the diffusivity of Si should be independent of (4H) × Si concentration. It would also, though, be the expected outcome for any simple hydroxylated vacancy-cation exchange such as the exchange of (2H) × Mg and Mg × Mg (Demouchy & Mackwell, 2003;Jollands et al., 2016), but would not be expected in situations where inter-site reactions are involved (Jollands et al., 2019). The exact pathway taken by the Si 4+ and 4H + in such an exchange is not clear, but it is likely that the Si 4+ moves temporarily into, or at least through, an interstitial site (Bejina et al., 1999;Brodholt & Refson, 2000).
Overall, the slow diffusion mechanism, slow , can probably be assigned to represent the diffusion of (4H) × Si , and thus is denoted Si herein. Note that this is an interpretation, but this is also supported by the agreement between the diffusion data from this study, and the Padrón-Navarta et al. (2014) data. However, it should also be noted that Padrón-Navarta et al. (2014) assigned their slowest diffusion mechanism to (4H) × Si based mainly on the observation that the 3,612 cm −1 band decreased. These (4H) × Si values are lower than those measured for H associated with M-site vacancies (Demouchy & Mackwell, 2003Jollands et al., 2016Jollands et al., , 2021Kohlstedt & Mackwell, 1998;Padrón-Navarta et al., 2014), which are themselves lower than Ds associated with H-2 H exchange (Du Frane & Tyburczy, 2012;Novella et al., 2017;Sun et al., 2019). The activation energy for the diffusion of hydroxylated M-site vacancies is lower (∼200-250 kJ mol −1 ), however, suggesting that at T > 1600°C, the hydroxylated Si-and M-vacancies may diffuse at the same rate. The (4H) × Si is 3-4 orders of magnitude higher than Mg diffusion in a H-free system, which is then 4-5 orders of magnitude higher than the equivalent Si diffusivity. This is generally lower than Si diffusion in an H-present system, although the magnitude of the effect of H is debated (see below). A compilation of these Ds is shown in Figure 9.

Assigning *
The process responsible for the elevated background absorbance of the 3,612 cm −1 band, described using the high diffusion coefficients, fast , is less straightforward. The H diffusion mechanisms that have diffusion coefficients high enough to homogenize (or nearly homogenize) a ∼1 mm olivine at the T-t conditions of these experiments are limited to the hydroxylated M-vacancy mechanism where (2H) × M exchanges places with Mg × M (Demouchy & Mackwell, 2003, proton-polaron exchange where H + moves by reducing Fe 3+ to Fe 2+ (Barth et al., 2019;Ferriss et al., 2018;Kohlstedt & Mackwell, 1998;Li et al., 2022;Mackwell & Kohlstedt, 1990) and H-2 H exchange (Du Frane & Tyburczy, 2012;Novella et al., 2017;Sun et al., 2019). Only the former is valid for pure forsterite in these experiments; the proton-polaron exchange, according to its current interpretation involving reduction of Fe 3+ , is irrelevant. Therefore, the elevated concentrations of the (4H) × Si defect associated with * fast probably forms by some interaction between mobile (2H) × M and some defect or species that is already present in the starting material.
Potentially, this could invoke trace elements, such as Al 3+ , notably present in the starting material at around 10-20 wt. ppm (see Zhukova et al. (2014) for analyses of the same starting material):

of 16
This is possible for these experiments, given that some bands associated with some Al-associated defects are present in some post-experiment FTIR spectra ( Figure 5). Another possibility could involve the removal of Si 4+ from its normal site into an interstitial position: This is an attractive possibility as it does not require any trace elements in the starting material, which may become important when considering Si diffusion experiments in extremely pure forsterite, discussed below. The Si •••• i defect is considered as possible in some early studies (Smyth & Stocker, 1975;Stocker & Smyth, 1978), but has recently been calculated to be more-or-less negligible (Muir et al., 2020), albeit in conditions of global equilibrium in a dry system, neither of which are appropriate here. Other reactions could be written but are unlikely, such as one which requires an energetically unfavorable Si vacancy initially: In the San Carlos olivine, there are more possibilities for forming (4H) × Si , given that there likely exists a faster "proton-polaron" exchange diffusion mechanism along with the M-site mechanism. Additionally, there likely exists Fe 3+ on the Si site (Nakamura & Schmalzried, 1983), which could be displaced to allow (4H) × Si to form. In summary, it seems reasonable to assume that the elevated (4H) × Si concentrations in the crystal cores relate to some form of diffusion-reaction process, where H is brought into the crystal by progressive (2H) × M ↔ Mg × M exchange, then interacts with some pre-existing defect to form (4H) × Si .

Notes Relating to Si Diffusion
An assumption underpinning the entirety of the following discussion is that, as stated in the introduction, silicon vacancies in H-free forsterite are unfavorable (Brodholt & Refson, 2000;Muir et al., 2022;Smyth & Stocker, 1975;Stocker & Smyth, 1978). It is assumed that whichever mechanism Si diffuses in the H-free system has a negligible contribution to Si diffusion in the H-bearing system. This is an end-member scenario based on ∼3 orders of magnitude difference in D Si between the H-free (slower) and H-bearing (faster) experiments of Dohmen et al. (2002) and Costa and Chakraborty (2008), respectively. However, comparing H-free D Si values from Fei et al. (2012) with H-bearing D Si from Fei et al. (2013) suggests a much smaller effect of the presence or absence of H, less than half an order of magnitude ( Figure 9). Whilst beyond the scope of this discussion, the intriguing effect of a thin ZrO 2 film on Si diffusion shown by Fei et al. (2012) (their Figure 6) might be worth revisiting. The presence of ZrO 2 leads to a considerable increase in D Si . The authors interpreted this as an experimental artifact, where a ZrO 2 film prevented shrinkage of an isotopically enriched Mg 2 SiO 4 layer, but another potential explanation is that the amount of ZrO 2 is sufficient to decrease the SiO 2 of the system to that associated with a Mg 2 SiO 4 -MgO-ZrO 2 assemblage, which might then considerably modify the defect population on the Si sublattice.
In any case, each time a (4H) × Si defect moves by (4H) × Si ↔ Si × Si exchange, as is interpreted to occur in the experiments described above, represented by slow (assigned to (4H) × Si ), a Si 4+ ion has to move in the opposite direction. The higher the concentration of the (4H) × Si defects and the higher their diffusivity, the greater the probability of (4H) × Si ↔Si × Si exchange in a given volume over a given time. Therefore, the higher the diffusivity of Si (D Si  (2017)). The measured absorbance was not converted to total absorbance (H. Fei, pers. comm. 2022), so this would need to be multiplied by ∼2 for this defect given the orientation of their samples (thick sections prepared || (010) Grant et al. (2007)). Additionally, they integrated over almost the entire O-H stretching region (3,000-4,000 cm −1 ), which appears to include a broad hump that should probably be removed by baseline subtraction prior to integration over a narrower range (e.g., 3,430-3,700 cm −1 ). Moreover, the Withers et al. (2012) absorption coefficients are likely more appropriate than those from Bell et al. (2003), given the spectra are dominated by the (4H) × Si defect, similar to the spectra in the samples analyzed by Withers et al. (2012). In fact, these corrections more-or-less cancel out, and their 810 wt. ppm H 2 O recalculates to around 940 wt. ppm. Regardless, their relationship between water content (denoted CH 2 O in their work) and D Si was D Si ∝ C 0.  Bell et al. (2003) calibration rather than the Withers et al. (2012), so again, these will cancel out to some extent. One important note is that their spectra from San Carlos olivine show strongest absorbance associated with a band at 3,573 cm −1 and a smaller band at 3,525 cm −1 , which can be attributed to a Ti-associated defect, Si } × , as described above. Resolving the spectra in their Figure 6 (middle spectra) into Gaussian components suggests that only around 30% of absorbance is associated with ( Because the Costa and Chakraborty (2008) and Fei et al. (2013) experiments were conducted at different pressures, the Si diffusion data cannot be directly compared. Whilst the activation volume for Si diffusion has been determined (Bejina et al., 1999;Fei et al., 2012), this was done for H-free olivine, and there is no reason to assume that the activation volume for Si diffusion containing (4H) × Si defects should be the same or even similar to the activation volume where these defects are absent. It is also possible to use the values of (4H) × Si from this study to predict D Si at a given . A water concentration of ∼5-10 wt. ppm H 2 O associated with (4H) × Si , as recalculated above for the 43 wt. ppm H 2 O experiment of Costa and Chakraborty (2008), is equivalent to on the order of ∼10 −4 to 10 −5 , also. At 1200°C, this study gives (4H) × Si of around 10 −13 , therefore log 10 (D Si [m 2 s −1 ]) is predicted to be ∼−17 to −18. This is ∼2 orders of magnitude higher than D Si measured by Costa and Chakraborty (2008), log 10 (D Si [m 2 s −1 ]) = −19.5 at 1200°C, 2 GPa. Whilst the data from this study suggest that the activation volume for (4H) × Si is minimal, it is probably inappropriate to extrapolate the data to 8 GPa, the pressure of the Fei et al. (2013) experiments, therefore a similar exercise is not done to predict D Si in their experiments from (4H) × Si in this study.
In both cases, the outcome is that the effect of (4H) × Si is smaller than would be predicted from the basic relationship between [(4H) × Si ], (4H) × Si , and D Si . In the first case, considering the 0.32 exponent determined by Fei et al. (2013), this mismatch is completely independent of the data from this study. This independence is important, because there still remains the (unlikely) possibility that there is a slower diffusion mechanism of (4H) × Si in olivine, which would not be seen in these experiments given the T-t conditions, that is, might be visible in longer/higher T experiments.

Explaining the Discrepancy With a Diffusion-Reaction Mechanism?
The explanation for elevated (4H) × Si concentrations in the crystal cores was that there exists some fast pathway for diffusion that forms (4H) × Si , without requiring any exchange between (4H) × Si and Si × Si . This concept may also help to explain why the D Si = (4H) × Si × [(4H) × Si ] relationship seems to fail here. In this relationship, it is implicitly assumed that every (4H) × Si defect moves, and specifically moves by (4H) × Si ↔ Si × Si exchange. If not, the defects are irrelevant to the calculation. Fei et al. (2013) show that for a ∼100× increase in [(4H) × Si ], D Si increases by ∼10×. This could then potentially be reconciled by proposing that only ∼10% of the additional (4H) × Si defects are moving by (4H) × Si ↔ Si × Si exchange at any given time. In order to reconcile the Costa and Chakraborty (2008) D Si with (4H) × Si from this study, we would similarly need to assume that only ∼1% of (4H) × Si defects move in this way.
Conceptually, it could be possible that rather than a (4H) × Si defect moving by exchanging positions with an adjacent Si × Si , the defect is destroyed, then reforms elsewhere, having moved via a different diffusion mechanism. This would lead to an apparently faster (4H) × Si mobility, and would have a considerably lesser effect on the mobility of To summarize, it is possible that some form of diffusion-reaction process is occurring. Whether this is sufficient to modify the effect of (4H) × Si on D Si to the extent described above is not clear. However, this does suggest that it is unlikely that Si diffusion can be described simply as (4H) × Si ↔ Si × Si exchange, even when (4H) × Si defects are present. This is consistent with the assignation of * fast to forming by some kind of diffusion-reaction process. It is also notable that Jollands et al. (2019) demonstrated that it is possible to remove almost all (4H) × Si defects from a metamorphic olivine at a rate which suggests that no/negligible (4H) × Si ↔ Si × Si exchange occurred.

Conclusions
The diffusion coefficients of H assumed to be associated with (4H) × Si ↔ Si × Si (i.e., (4H) × Si ) exchange have been determined in pure forsterite and natural olivine at 1200°C-1400°C. These are in good agreement with a previous determination of diffusion coefficients assigned to the same mechanism. These results confirm that the diffusivity of the hydroxylated Si vacancy is lower than that of the M-site vacancy, with a higher activation energy. Attempts to link the determined (4H) × Si to Si determined in H-bearing systems were unsuccessful, suggesting that simple assumptions of Si diffusion by (4H) × Si ↔ Si × Si exchange are likely inappropriate.