The effect of water and pressure on fabric development in olivine

Water has a very strong effect on both the strength and fabric development of forsterite but the mechanism of this effect is unclear. In the paper we use Density Functional Theory Peierls-Nabarro modelling to examine the effect of water on the Peierls stress of different forsterite slip systems. We find that water in Mg vacancies will weaken [100](010) slip and thus produce A fabrics while water in Si vacancies weakens [001](100) slip and thus produces C fabrics. With a combination of DFT and forcefield based methods we find that while in the bulk hydrated vacancies typically occur as M1>M2>Si, near the [100] and [001] screw dislocation cores they occur as Si>M2>M1. Thus by simple modification of the Peierls stress and by segregation of Si vacancies to dislocations, water enhances C fabrics in forsterite. This production is enhanced by pressure and should be independent of water concentration above ~90 ppm H/Si in perfect forsterite. The shape of dislocations is not significantly modified by water and thus the stress exponent of dislocation glide should be unaffected. Water additionally reduces the activation volume of the Peierls mechanism by ~5 times by changing the pressure dependence of slip. E fabrics cannot be produced by a dislocation creep mechanism due to the consistently high Peierls stress of the [100](001) slip system and thus require another deformation method to form.

significantly modified by water and thus the stress exponent of dislocation glide should be unaffected. Water additionally reduces the activation volume of the Peierls mechanism by ~5 times by changing the pressure dependence of slip. E fabrics cannot be produced by a dislocation creep mechanism due to the consistently high Peierls stress of the [100](001) slip system and thus require another deformation method to form.

Introduction:
Olivine is a major component (40-80%) of the Earth's upper mantle and thus understanding its deformation is critical for understanding upper mantle rheology. The plastic deformation of olivine, however, is very complex and not fully understood and the effect of water on this deformation adds further complexity. Olivine can dissolve large quantities of water (up to ~30,000 ppm H/Si depending on pressure (Hirschmann, 2006)) with the water content of the deep asthenosphere estimated from conductivity data to be around 1,500±800 ppm H/Si (Wang et al., 2006, Karato, 2011 but with heterogeneous distribution of between 0-10,000 ppm H/Si (Yoshino et al., 2006, Karato, 2011.  (Raterron et al., 2009, Raterron et al., 2007, Raterron et al., 2011, Raterron et al., 2012, Bollinger et al., 2016, Jung et al., 2009, Ohuchi et al., 2011. Similarly the introduction of water has been shown to induce the formation of E or C fabrics (Jung and Karato, 2001, Jung et al., 2006, Katayama et al., 2004, Ohuchi and Irifune, 2014. All of these papers show, however, that fabric transitions in olivine possess a complex interplay between water content, pressure and stress/temperature that is not yet fully understood. Concluding information about deformation mechanics from natural rocks is difficult due to a lack of information about their deformation environments and histories but in general there is a weak trend of C fabrics being found in wet rocks (which supports water inducing an A to C fabric change) and B fabrics being found in rocks from convergent boundaries (which supports stress inducing an A to B fabric change) while rocks with E fabrics are relatively rare (Karato et al., 2008b, Michibayashi et al., 2016, Jung, 2017. Understanding texture development in olivine is important because different olivine fabrics have different seismic anisotropies and thus changes in olivine fabric preference will change the flow properties of the upper mantle. The presence and the change of olivine fabrics has been used to explain various properties of the upper mantle including the decrease in seismic anisotropy and the sign change of VSH-VSV with depth (both indicative of an A to B/C transition) (Ohuchi and Irifune, 2013, Ohuchi and Irifune, 2014, Raterron et al., 2012, Ohuchi et al., 2011, the larger anisotropy in continental vs oceanic lithosphere (related to the depth of the A to B transition) (Raterron et al., 2012), the conversion of shear wave splitting from trench parallel near the trench to trench perpendicular near the back arc in colder wetter slabs (conversion from A/C to B) (Jung et al., 2006, Karato et al., 2008a, Ohuchi et al., 2012, Katayama and Karato, 2006, Kneller et al., 2005, Kneller et al., 2007 and finally the overall anisotropy and shear wave splitting direction of the upper mantle (predominately A or E fabric) (Karato et al., 2008b). Thus knowledge of the effect of water on these transitions is essential to constraining the amount of water in the upper mantle, its overall distribution and the anisotropy and thus the flow of the upper mantle.
The mechanism for how water changes fabric texture in forsterite, however, is unclear as is the interaction of this mechanism with pressure and stress. A major complication in determining this mechanism is the complex nature of forsterite deformation with upper mantle conditions being conducive to both dislocation and diffusion creep with the latter favoured by lower temperatures, smaller grains sizes and lower stresses (Faul et al., 2011, Ohuchi et al., 2017, Fei et al., 2016, Katayama and Karato, 2008, Wang, 2016, Kawazoe et al., 2009). Dislocation creep is likely dominant in the asthenosphere and diffusion creep in the lithosphere (Fei et al., 2016) but both deformation mechanisms will likely carry strain under most conditions (Nishihara et al., 2014) and dislocation glide is possible at even quite low temperatures and stresses (~800 K <200 MPa) (Boioli et al., 2015b). In this paper we shall address the viability of one of the most straightforward ways that water could affect fabric evolution-by changing the Peierls stress (ease of dislocation glide) of different slip systems (Mackwell et al., 1985, Katayama andKarato, 2008) and thus changing which slip systems carry the most strain. We shall do this by using atomistic Peierls-Nabarro calculations to model the different slip systems and to calculate their Peierls stress. Previous work (Skelton and Walker 2017)  and [001](100) slip systems while also considering the effect of pressure which is experimentally hard to constrain while also controlling water fugacity.

Methods
To establish the effect of water on fabric development in forsterite two things need to be established-the position of the water in the crystal and its effect on the forsterite slip systems. To establish the former we calculated the energies of dry and hydrated Mg and Si vacancies in a forsterite bulk (using Density Functional Theory (DFT)) and near a dislocation core (using forcefields). To examine the latter point we calculated the Peierls stress of To calculate the energies of vacancies and the elastic constants a (2x1x2) forsterite super cell was used with a plane wave cut off of 1000 eV, (4x4x4) k-points and relaxed to a force tolerance of 0.01 eV/Å and an energy tolerance of 1x10 -5 eV/atom. A (2x1x2) supercell was used to ensure that there was roughly 10 Å between repeating vacancies in all directions.
A (4x2x4) supercell was also tested and the vacancy energy changed by <0.01 eV. While hydrated vacancies are chargeless, dry vacancies are charged and so the energy produced by DFT includes a defect-defect interaction term while we require the energy of a defect in infinite space which would not include such a term. We can correct for this by assuming our system is an array of point charges in a neutralising background charge as outlined in Leslie and Gillan (1985) and previously used for forsterite by Brodholt and Refson (2000). To use this method the relative permittivity of the cell needs to be set -we used a value of 6.2 following Brodholt and Refson (2000).

Vacancy Segregation
To calculate the energy of vacancies at a dislocation core it is necessary to model the core which requires 10,000s of atoms and thus is beyond the scope of DFT. Instead we use interatomic potentials to calculate the energy of the segregation of vacancies from a bulk like environment (where their energy is calculated by DFT) to the region of a dislocation core. To do this we used the same methods and core structures as work (Skelton and Walker 2017) but we modelled Si vacancies rather than Mg vacancies. This method is outlined in more detail in the supplementary methods and in work (Skelton and Walker 2017) but in essence a dislocation core was constructed using the disloPy code (https://github.com/andreww/disloPy/) and then the energy of a dry or hydrated Si vacancy is calculated at all possible sites relative to the dislocation core up to a set distance away (12.5 Å). Energies were calculated using a forcefield designed for hydrated olivine , Lewis and Catlow, 1985, Schroder et al., 1992 (Table S1) and the GULP code (Gale, 1997).

Peierls-Nabarro Method:
To calculate the dislocation structures and Peierls stress we use the Peierls-Nabarro method . This is an approximate method which represents a dislocation as a series of partial dislocations across a glide plane with the assumption that the dislocation is collinear. The dislocation is then composed of two forces-the linear elastic response of the crystal and the energy of inelastic misfit-which balance at equilibrium. The traditional PN formulation requires a continuous dislocation density which is computationally difficult to calculate and so instead we use a semi-discrete PN formulation  which allows us to calculate the dislocation density at a series of isolated points which can be run as isolated computations. The PN calculations were done using the disloPy code (https://github.com/andreww/disloPy/) which has been previously documented in Skelton and Walker (2017)  An additional consideration concerns the effect of breaking the symmetry of the hydrated vacancy during displacement and thus introducing entropy barriers. This concern is addressed in the supplementary information and is found to be negligible.

D-Rex Calculations
To produce illustrative fabric diagrams we used the D-Rex approach (Kaminski et al., 2004) as implemented by C. Thissen (https://github.com/cthissen/Drex-MATLAB/). This is a kinematic model for simulating texture development through considering energy reductions of an averaged field (rather than individual grains) using dynamic recrystallization, subgrain rotation and grain-boundary migration which are all set parameters. In all cases 1000 particles were used, final shear was set to 1 and our velocity gradient tensor was proper exploration of these values will be required. It is also required to set the critical relative shear stress (CRSS) for the four different slip systems. For each of our calculations these were set to the screw Peierls stresses at the appropriate conditions that we determined through our calculations. The stress exponent was fixed to 3.5 for all slip systems.

Hydrated Vacancies in Forsterite:
Before we can comment on the effect of water on forsterite slip systems we need to know where the water exists in the crystal. Water distribution in olivine and forsterite is a complex problem but 4 major sites of water have been determined that are (in Kroger-Vink notation) (2 ) , (4 ) , (1 ) ′ · and ·· (2 ) ′′ (see for example Tollan et al. and Si (4 ) vacancies. There are a few reasons for this choice. First (2 ) and (4 ) vacancies will always be present in forsterite whereas (1 ) ′ · and ·· (2 ) ′′ require contaminants (and in the case of (1 ) ′ · high oxygen fugacity) to form and thus we can establish the baseline effect of water on perfect forsterite. Second at high water concentrations which promote dislocation glide (Katayama and Karato, 2008) water is likely to outnumber the concentration of Ti and Ferric iron defects. Third we calculated the binding enthalpy (the enthalpy to remove the hydrous component from the non-hydrous component) of these two latter defects to be high-2.2 and 5.0 eV for (1 ) ′ · and ·· (2 ) ′′ respectively at 5 GPa and static conditions. This means that the hydrous vacancy part of these complexes is unlikely to segregate to an forsterite dislocation (as the segregation energy of hydrous vacancies shown below is generally lower than these binding energies) and thus the entire complex would need to be segregated to the dislocation. Thus while hydrated vacancies segregate easily to dislocations as discussed below (1 ) ′ · and ·· (2 ) ′′ may not segregate so easily though this needs further study. For this study the only important hydrated complexes are those that will segregate to the dislocation core. Thus we shall examine (2 ) and (4 ) (and ′′ and ′′′′ ) which shall allow us to examine the dislocation mechanisms of perfect forsterite, very wet forsterite and likely also mildly wet forsterite.
Therefore we considered water existing as hydrated Mg vacancies (2H + + ′′ ) on the M1 or the M2 site or as hydrated Si vacancies (4H + + ′′′′ ) and calculated the relative energy for each. The relative energy of M1 and M2 hydrated vacancies is simply their energy difference, the relative energy of M1 or M2 vacancies to Si vacancies depends upon the background chemical environment. There are two major cases: 2(2 ) 4 + 2 2 4 → 2 (4 ) 4 + 2 2 2 6 Reaction 1 2(2 ) 4 + 4 → 2 (4 ) 4 + 2 2 4 Reaction 2 Where the first reaction represents a perioditic upper mantle where orthoenstatite is stable and reaction 2 represents a silica depleted composition. The same reaction occurs in the dry case but with dry vacancies.
These two reactions have been previously considered in detail using theoretical methods (Walker et al., 2007, Qin et al., 2018 where it is found that in enstatite (2 ) vacancies are largely favoured but in periclase (2 ) and (4 ) have similar concentrations with exact conditions being important. We calculated the energies of Reaction 1 and 2 (Table 1) and found largely similar answers to these previous works with enstatite distributions shown in Figure 1 and periclase distributions in Figure S1. While (4 ) is energetically favoured over (2 ) the large configurational entropy gains of pushing Reaction 1 and 2 to the right see (2 ) vacancies dominate in enstatite-buffered conditions in the lower mantle with both increased temperature or pressure favouring (2 ) and increased water concentration favouring (4 ) . In the case of dry vacancies Mg vacancies are overwhelmingly favoured in all conditions as both the energy and configurational terms strongly favour their formation (Table 1). In both the wet and dry case M1 Mg vacancies are significantly favoured (>0.6 eV) over M2 Mg vacancies (Table 1).
In the case of dislocation glide, however, what is important is not vacancies in a bulk crystal but vacancies that are near a dislocation which is a very different chemical environment. Simulating very large dislocations with DFT is unfeasible and so instead forcefield methods are required. The segregation of Mg vacancies in forsterite clusters was previously published in work (Skelton and Walker 2017) and here we extend this method to Si vacancies. Figure S2 shows the segregation of dry Si vacancies and Table 2 lists the maximum energy gain of moving dry and wet vacancies to dislocation core taking into account the configurational entropy that is lost by such a movement. The energy of hydrated Silicon vacancies that is calculated by forcefields is extremely sensitive to the starting position of the hydrogen. Thus rather than generate segregation maps like we did for dry silicon instead we focused on getting the most accurate segregation energy for hydrated Si by sampling a vacancy removed from the dislocation and the 3 vacancies closest to the dislocation core and in each case trying 6 different starting hydrogen arrangements. This then allows us to calculate a segregation energy which is reported in Table 2. While the absolute energy of Si vacancies and particularly hydrated Si vacancies is likely poorly represented by forcefield calculations the relative energy of changing the stress field in which the vacancy sits should be much better formulated and thus the segregation energy of the vacancies should be reasonably accurate.
In all cases there is a strong enthalpy gain in moving a vacancy to a dislocation core with M1 vacancies gaining ~-1 eV and M2 and Si vacancies ~-2-3 eV and a strong entropy loss which can vary from +0.5-4+ eV depending on temperature and concentration but this loss is much lower for Si vacancies. As shown in Table 2 these effects combine such that Si vacancies strongly segregate to dislocation cores whereas M1 vacancies remain in the bulk and the position of M2 vacancies is condition dependant. At the dislocation cores, however, hydrated Si vacancies will be strongly favoured. Pressure, temperature and water content have little effect on this partitioning such that at all upper mantle conditions hydrated Si vacancies should dominate at dislocation cores. Figure 1 and S1 show that at decreasing temperature, decreasing SiO2 activity and increasing water content the dominant water vacancy becomes hydrated Si vacancies at dislocation cores. Due to this we shall focus on the effects of hydrated Si vacancies on the Peierls stress but we shall also discuss the effects of any M1 or M2 vacancies that are present at dislocation cores.

Dislocation Glide in Forsterite
Elastic Constants: Our calculated athermal elastic constants are presented in Table 3 alongside the anisotropic energy factor for the four different slip systems. These are similar to those calculated by Durinck et al. (2005). Elastic energies of a perfect crystal were used for all systems including those with vacancies as they more reliably reflect the elastic forces between individual atoms and as hydration has little effect on the elasticity of forsterite (Cline et al., 2018).

Base Forsterite
In order to consider the effect of water we first need to consider the deformation of a crystal with no vacancies. The slip planes which produce the lowest energy GSF profiles (as determined by manual search) are shown in Figure 2 -in all cases slip planes are more favourable if they avoid severing Si tetrahedra though this is not possible in the [100](001) slip system. Sample GSF and dislocation profiles can be seen in Figure 3 and 0 GPa Peierls stress and dislocation widths in There is some controversy over the [100](010) slip system as Mahendran et al. (2017) found it has a tendency to dislocate into two partial dislocations but this was not observed by us or by Durinck et al. (2005). Dissociation of a dislocation core is typically indicative of a Peierls stresses) because its slip plane passes through the centre of a SiO4 tetrahedra which is highly unfavourable.
As shown in Table 4 (Raterron et al., 2009, Raterron et al., 2007, Raterron et al., 2011, Raterron et al., 2012, Ohuchi et al., 2012. This discrepancy is due to the non-planar nature of the [001] screw dislocation (Mahendran et al., 2017, Carrez et al., 2008. Using interatomic forcefields on large systems and the application of direct shear to the system Mahendran et al. (2017) found that in order to glide [001] screw dislocations must first convert to a planar core. The results of Mahendran et al. (2017) are shown in Table 4 and compared to our results. There are two significant differences between our computational regime and that of Mahendran et al. (2017)-one in computational approach (DFT vs forcefield) and one in method (PN vs directly applied stress). Our approach should determine the resistance to flow of a planar core more accurately but it does not consider the rate of conversion of non-planar cores to planar cores and so we need to apply a correction to our  (010) and [001](010) slip systems are less affected. This is because in the former two systems the slip planes cuts through the centre or the corner of a SiO4 tetrahedral and so replacing the Si with a vacancy reduces the barrier to distorting/breaking the tetrahedron. The difference between hydrated and non-hydrated vacancies (shown in Fig S10) is small for the most stabilised [001](100) system but is on the same order of adding a vacancy for the other 3 slip systems. Dry Si vacancies are exceedingly rare (Table 1), however, and so only wet Si vacancies need to be considered.  Figure S11-S12).
[100](010) slip is easiest at all pressures but the barrier to planar glide of the [001](010) slip system approaches 0 when M2 vacancies are present and its unfavourability is entirely due to its assigned locking-unlocking rate. For those systems with (010) slip planes the large reduction in Peierls stress with M2 vacancies is due to the slip planes containing an M2 atom and removing this atom removes a barrier to slip which occurs similarly in both wet and dry cases. In the other two slip systems there are complex interplays between charge, hydrogen positioning and slip plane barriers.
M1 vacancies are unlikely to be present at dislocations but we shall still consider their effects ( Figure S13 and S14  (010) plane and so wet and dry M1 vacancies operate similarly.

Fabric Development in the presence of Vacancies
To illustrate these cases we then predicted the resulting fabric of our forsterite in various conditions using the D-Rex code. Using our raw uncorrected data in all tested cases a clear B fabric develops ( Figure S15). This is expected as the Peierls stress of the [001](010) slip system is the lowest when not considering the effects of [001] non-planarity. With our corrected data (Figure 6 and 7) we see that water can make clear modifications to the texture.
In the absence of vacancies, increasing the pressure switches the developed fabric from A to B. With hydrated Si vacancies (non-hydrated vacancies are shown in Figure S16) an A texture with some weak C texture is seen at 0 GPa that develops by 10 GPa into a mixture of A, B and C textures. This effect happens because while the [001](100) and [001](010) slip systems have Peierls stresses that mildly harden/weaken with pressure ( Figure 4) the [100](010) slip system hardens significantly with pressure and thus the A fabric weakens with pressure. This indicates that some pressure, as well as water, is needed to form C fabrics in forsterite. With hydrated M1 and M2 vacancies (Figure 7) a strong A fabric is observed throughout our pressure range-as expected from their slip system activities-and so only Si vacancies can induce strong texture changes in forsterite.

Discussion:
Our main overall conclusion is that water induces a C fabric and somewhat a B fabric In dry (<40 ppm) forsterite a fabric transition from A to B (or B-like with some C components) occurs with pressure. This has been observed with EBSD measurements after deformation (Jung et al., 2009, Ohuchi et al., 2011, Bollinger et al., 2016 and via direct monitoring of the strength of different slip systems (Raterron et al., 2009, Raterron et al., 2007, Raterron et al., 2011, Raterron et al., 2012. The pressure of this transition is varied from ~3 GPa (Jung et al., 2009) to 6-8 GPa (Ohuchi et al., 2011, Raterron et al., 2007, Raterron et al., 2011. We find that this transition can be easily explained by the different pressure derivatives of the slip systems as has also been seen by Raterron et al. (2011 screw core. This rate (which we cannot determine with these methods) is highly likely to be controlled by stress. A stress-induced decrease in the locking-unlocking rate would explain the stress induced transition to B fabrics that has been previously observed (Jung and Karato, 2001, Jung et al., 2006, Katayama and Karato, 2006 at ~300 MPa. Additionally in Ohuchi et al. (2011) it was observed that B fabrics (as opposed to mixed B-like fabrics) formed with increasing clarity at lower temperatures but these lower temperature runs also had higher stresses which would aid B fabric formation under such a system.
In wet forsterite a variety of effects have been seen. Firstly at low stresses <~200 MPa and pressures ~2 GPa an A to E transition was observed for moderate amounts of water (< 800 ppm H/Si) and an A to C or E to C transition for larger amounts of water (Jung et al., 2006, Katayama et al., 2004, Jung and Karato, 2001. At higher pressures (>10 GPa) a C fabric is developed with a very wide range of water (50-29000 ppm H/Si) (Couvy et al., 2004, Ohuchi andIrifune, 2014). Our results predict that water will promote C (and somewhat B) fabrics simply by hydrous modification of relative Peierls stresses but that some additional pressure is required to resolve the fabrics which explains the fact that C fabrics are more commonly observed at high pressures in experiments.
The production of E fabrics with moderate water contents is more vexing. Our results find that under all pressure and vacancy conditions the [100](001) slip system has a considerably higher Peierls stress than other slip systems. There is no obvious modification that will lower the Peierls stress of this slip system as it will always sever a SiO4 tetrahedron which is energetically unfavourable. Experimentally Tielke et al. (2016)  is ~200 ppm (H/Si) for dislocations with (010) slip and ~300 ppm for those with (100) slip. The required concentration to induce a fabric is likely much lower than this as ~10% of the dislocations carries 95% of the strain in forsterite (Boioli et al., 2015a) which would move this transition to ~60-90 ppm. Beyond this saturation point the water will overwhelmingly reside in the bulk and should have no effect on the fabric transition and thus for any reasonably wet (>~90 ppm) sample water should simply enable C fabrics with no concentration effects. Under more high stress conditions (450 MPa) we find that that 10% saturation would be achieved at a water concentration of 200-300 ppm and thus stress has a small effect on the saturation limit of dislocations.
A different set of relations has been seen with wet forsterite at high stress. At low pressures (~2 GPa) stress above ~300 MPa has been seen to convert A/E/C fabrics to B fabrics (Katayama and Karato, 2008, Jung et al., 2006, Jung and Karato, 2001, Katayama and Karato, 2006, and at high pressures (7.2-11.1 GPa, 400 MPa stress) (Ohuchi and Irifune, 2013) to convert C fabrics to B and then to A fabrics with increasing amounts of water though there is some debate about this barrier as B fabrics have been produced in wet forsterite at ~2.1-5.2 GPa and below <300 MPa stresses (Ohuchi et al., 2012). As temperature decreases the required stress to convert to B fabrics will also decrease (Katayama and Karato, 2006)-at a typical lower mantle stress of 80 MPa C fabrics will be produced in favour of B at temperatures greater than 1173 K but at a high stress of 450 MPa this temperature rises to ~1600 K. As discussed above the conversion of A or C fabrics to B fabrics is likely related to the unlockinglocking rate of the [001] screw dislocation. In both the dry and the wet case planar glide of the [001](010) slip system is the easiest slip system and so increasing the rate at which this screw dislocation becomes planar (if the locking-unlocking rate decreases with stress) will improve the formation of B fabrics at all water concentrations. The formation of A fabrics with large amounts of water and stress (Ohuchi and Irifune, 2013) has a less obvious source but is perhaps related to the fact that Mg vacancies (which will be present with such large amounts of water) promote A fabrics.
There are two other common stress mechanisms but they are unlikely to cause the formation of B fabrics. The most obvious effect of stress on dislocation creep is an increase in the number of dislocations with dislocation density increasing from 5.3E+11 m -2 at 80 MPa to 6.0E+12 m -2 at 450 MPa (Karato and Jung, 2003). The variation in these values is much lower than the variation of water contents across experimental ranges which does not induce B fabrics and so this should not be a factor in B fabric formation. The other possibility is that stress has a varying effects on the ease of glide of different slip systems. This will not be a factor if the stress exponents are the same for all slip systems which is true for dry forsterite (Durham et al., 1977) but is possibly not true for wet forsterite. We find that water does not have a large effect on the shape of dislocations and thus it is unlikely to affect the stress exponent of the dislocation and thus we expect wet forsterite to have the same stress exponent for all slip systems.  (010), (100) and (110) planes as well as (001)[010] slip (Demouchy et al., 2013, Durham et al., 1977, Phakey et al., 1971, Demouchy et al., 2009, Evans and Goetze, 1979, Raleigh, 1968). Polycrystalline powders also show a rapid conversion from low temperature to high temperature dislocation glide patterns at ~1373 K (Demouchy et al., 2013, Bai et al., 1991, Raleigh, 1968 (Bai et al., 1991) and thus more work would be needed to evaluate these systems which would include considering diffusional effects and calculating the slip of non-standard glide planes.
A final consideration is the effect of water and pressure on the activation volume of forsterite which is important for calculating its deformation as a function of depth. In dry olivine the activation volume has been observed to decrease from 15±3 to 0±1.2 cm/mol 3 by increasing the pressure (Raterron et al., 2011) Korenaga and Karato, 2008) found that dry olivine has an activation volume of 13±8 and wet olivine 4±3 such that a significant decrease in activation volume with hydration is likely. The large scatter of activation volumes of nominally dry olivine is likely also an indicator that activation volume is sensitive to water fugacity (Ohuchi et al., 2017). We find that the pressure dependence of the Peierls stress decreases by ~5 times (depending on pressure) (         Table 1: Energy (in eV) of converting an M1 vacancy to an M2 vacancy or Reaction 1 and 2 in dry or wet conditions with positive values representing the formation of M1 vacancies. The first three columns show the variation in enthalpy of the reaction as a function of pressure and the last three show the energy of the reaction at a specific temperature and pressure by adding the enthalpy to the calculated configurational entropy. These configurational entropy effects are invariant to pressure and scale with temperature. Vacancy concentrations are given in ppm H/Si-for dry vacancies which contain no H they are set to equivalent concentration of the wet vacancies.    Table 5: Peierls stress of forsterite in perfect and wet (containing a hydrated Si vacancy) conditions and the pressure derivative of this stress as determined by fitting a polynomial to these points.
Supplementary Method: Peierls-Nabarro: To calculate dislocation glide in forsterite we use the Peierls-Nabbaro (PN) formalism . In this model a dislocation is represented as a distribution of partial dislocations along the glide plane. Two different forces determine the shape of the dislocation: the elastic interaction energies between atoms (which broaden the dislocation) and the inelastic misfit energy caused by the introduction of disregistry by forming a dislocation (which narrows the dislocation). At equilibrium these two forces balance and a dislocation width can be determined-this is the key output of the PN model. Once this has been established the Peierls stress can be determined.
Our Peierls-Nabarro (PN) calculations were done using the disloPy code (www.github.com/andreww/dislopy). This formulation represents a planar dislocation with finite core-width as a distribution of dislocation density ρ along the glide plane (see Bulatov and Cai (2006) for a complete treatment). The classical PN treatment requires a continuous dislocation density and misfit energy. This has two issues-1) it is translationally invariant and so mobile under an infinitesimal external stress (and thus the Peierls stress is impossible to calculate) and 2) it is computationally difficult to calculate. Instead we shall use a semidiscrete formation  where the dislocation is represented by a distribution of ρ on a discrete lattice and the inelastic energy is computed by summing over this lattice. This has the distinct advantage of being able to calculate the inelastic energy at a series of individual lattice points which can be done as a series of independent atomistic calculations rather than having to calculate the entire dislocation in a single calculation.
There are three key assumptions in this model-firstly that the non-linear interactions between adjacent partial dislocations are negligible, secondly that the core is relatively compact and thirdly that it is collinear. While the first two are likely to be true in forsterite the third is not for [001] screw dislocations-the breakdown of this assumption shall be revisited in the results and discussion sections.
To represent the two different forces two different functions need to be calculated in two separate calculations. To calculate the elastic interaction we need K, the elastic prefactor, which can be determined from the elastic constants of the appropriate crystal. These are calculated for forsterite with Density Functional Theory (DFT) as explained in the text. As forsterite is inherently anisotropic we then use the sextic formulation of Stroh (1958) to convert the elastic constants into an appropriate K for each slip system.
To determine the misfit energy we need to know the energy value of the function γ(u). This is called the γ-line in one direction and the γ-surface in two directions and is also called the Generalised Stacking Fault (GSF). This function gives the inelastic misfit energy of a crystal along the discrete lattice γ offset by the disregistry (u-the displacement of atoms from their perfect crystal positions). This is calculated with γ being set to the required dislocation vector in the required dislocation plane. To calculate the GSF we use a different set of DFT calculations. To do this we take a slab with long repetitions normal to the shear plane and then displace half of it various directions along γ and allow atoms to only relax normal to the shear plane. For each distance along γ the energy of the displaced slab is obtained, the energy of the crystal without displacement is subtracted, the remaining energy is divided by the surface area of the plane and then a γ(u) function is constructed.
With these two functions calculated the PN calculation proceeds as follows.
Firstly the dislocation density (ρ) is written in terms of the disregistry (u): Where x is the coordinate along the gamma surface (displacement direction of the dislocation).
The energy of the dislocation is then represented with Equation 2: Where Emisfit is the inelastic energy, Eelastic is the elastic energy and Ework is the work done on a dislocation by an applied stress (σ). These terms are then represented by: = ∑ ( ( )) Equation 5 In equation 4 K is the elastic prefactor defined above. In equation 5 ap is the spacing between adjacent atomic planes and n are the grid points where calculations were done.
By setting Equation 3 to 0 (σ=0 no applied stress) the equilibrium disregistry function u(x) (static core structure) can be calculated by minimising Ecore in Equation 2. This minimisation is constrained by the fact that the sum of ρ(x) must equal the burgers vector b. To find the Peierls stress you then progressively increase σ and relax the disregistry at each step, eventually no minimum will be able to be found-ie the dislocation is completely mobileand this will give you a Peierls stress (σ).
To solve equation 2 numerically the disregistry function u(x) needs to be found. We do this by representing it as the sum of arctargent functions (ie partial dislocations) using To calculate the segregation of vacancies to clusters we first need to model a dislocation. This is too large to be modelled with DFT and so instead we use a forcefield based method. All segregation calculations were done with GULP (Gale, 1997) and the a forcefield designed for wet olivine . Construction of the dislocation cores and of the segregation maps was done with the disloPy code (www.github.com/andreww/dislopy).
The potential that we used models cations with a formal charge (Mg 2+, Si 4+) whereas O atoms are modelled as a positively charged core (+0.84819) with a negatively charged massless shell (-2.84819). All cation-anion pairs are joined by a Buckingham potential and SiO4 tetrahedra are fixed with a harmonic three body potential with details of this given in Table S1.
To prevent dislocation-dislocation interactions which are often large we instead use a 1D periodic cylinder which is periodic along the dislocation line and contains the dislocation in its centre. This cylinder has two radii-the first (R1) contains the dislocation and its surrounding atoms, the second (R2) contains atoms fixed in their bulk positions to simulate a perfect crystal.
Dislocation cores were constructed as in work (Skelton and Walker 2017)with the same cores and dislocation centres found. Dislocation core energies were found to converge (within 10 meV/ Å 2 ) at R1= (Skelton and Walker 2017)-their most stable cores were found to be at [0.74775,0.5] and [0.25,0.50225] and the energy of their cores converged when R1 was set to 25 and 35 Å respectively. Columbic energy was calculated using the Wolf sum (Wolf et al., 1999) which generally gives more physically reasonable properties for systems without 3D periodic boundaries and its cutoff was set to 15 Å with a damping parameter x=0.2 Å -1 . The effect of the atoms fixed in their bulk positions in the outer circle cannot extend beyond the columbic cutoff and so R2 was simply set to R1+15 Å. The cluster needs to be repeated along the periodic dimension and this was done 4 times for Si vacancies with dislocation line vector [001] and 5 times for those with dislocation line vector [100]. The minimum distance between Si vacancies and their periodic repeat unit was 23.91 Å and 23.95 Å in the [100] and [001] orientated dislocations respectively.
The segregation energy of a defect is defined as the energy of reaction from taking it from a site in the bulk of the crystal and placing it near a dislocation core. In other words: Where ΔEdefect is the excess energy of the defect in the bulk and ΔEsite-bulk is the excess energy of a site compared to a site in the bulk. ΔEsite-bulk is calculated simply by removing a Mg or Si from each possible site in a different run and then comparing the energy of the different runs and ΔEdefect is calculated from 3D periodic boundary calculations.

Configurational Entropy of Vacancies
There are two major forms of configurational entropy in olivine. The first is site dependant and simply depends on the ratio of Mg vacancies to Mg sites or Si vacancies to Si sites. The second only occurs in hydrated olivine and concerns the arrangements of hydrogen in the vacancy.
For the distribution of Si vacancies in the bulk we can simply use Boltzmann's entropy: Where W is the amount of ways the system can be arranged. To solve this for Si vacancies in the bulk we used Stirling's approxmation Where N is the total number of Si sites and nvac is the number of Si vacancies.
For Mg vacancies there are two non-equivalent sites and thus we must use the Gibbs entropy formulation: Where pi is the probability of the occurrence of each microstate i. This is determined by Where Z is the canonical partition function

= ∑ (− ) Equation 12
For the different vacancy sites the energy of the microstates with vacancies on M1 or M2 respectively were taken from Table 1. For vacancies at dislocation cores they are pinned to a specific site and so this term was set to 0. Strictly there is a multitude of sites near dislocation core and their energy should be formulated, their effect on dislocation glide tallied and their distribution calculated but this is a simplified model that captures the partitioning of vacancies to the most stable dislocation core site and the effect of this vacancy on the core glide dynamics.
To calculate the configurational entropy of hydrogen in the vacancy we again use a Gibb's distribution (Equation 10-12). The microstates and their energies in this case were taken directly from Qin et al. (2018).
The effect of configurational entropy on Peierls-Nabarro calculations: An additional consideration when dealing with hydrated vacancies is the loss of symmetry and thus configurational entropy caused by displacements (and in some case destruction) of the MgO6 and SiO4 tetrahedra. In the general sense this is dealt with by letting hydrogen find the minimum energy arrangement. As the Peierls-Nabarro model does not have a temperature component this should be sufficient. Real displacements happen at temperature, however, and breaking the equivalence of various hydrogen arrangements may incur an energy penalty, through loss of configurational entropy, not adequately accounted for by the model.
To examine the effect of this we took a sample case-[100](010) slip with an M2 vacancy. This was chosen because the M2 octahedra containing the vacancy is split during displacement and thus this system will show the maximum effect of configurational entropy. With this case we applied our PN model in two different ways. Firstly we did a standard GSF calculation which assumes no change in configurational entropy throughout and allows the H atoms to relax freely. Secondly we calculated the γ-line maximum and minimums of all 9 hydrogen arrangements (with hydrogen fixed to a specific O site but allowed to relax in that site) and then created a single γ-line maximum and minimum by using the Gibbs entropy formula (Equation 10-12) to partition the energy of all 9 individual pathways.
What was found was that displacing the crystal causes the most stable H configuration to be considerably more favoured than when compared to the undisplaced cell. At the γ-line minimum (undisplaced unit cell) the energies from different H configurations range from 0 to 0.83 eV with the nearest configuration to the minimum energy configuration being at 0.05 eV-nearly identical. At the gamma-line maximum (where the unit cell is at its most distorted) the range of energies for different H configurations now range from 0-4.48 eV with the nearest energy to the minimum configuration being at 0.68 eV. This vast discrepancy is because all of the H configurations except the most stable have H bonds across the shear plane and thus are strongly disrupted by shear. This means that when considering the statistical distribution of H arrangements, the most stable arrangement is overwhelmingly dominant and this is what is calculated in the regular GSF calculation. The difference in in GSF maximum from our standard calculation and from a calculation with explicit entropy consideration was found to be less than 5 meV/ Å 2 at 1500 K which is about equal to the error we obtain from not having infinite layers. Thus this is not a significant effect and shall be ignored. Figure S1: Ratio of the concentration of wet Si vacancies compared to wet Mg vacancies (M1+M2) for a varying water content at 5 GPa, in the presence of periclase (Reaction 2) and at different temperatures (1000, 1500 and 2000 K blue green and red respectively). Three different cases are shown-one with a pure bulk system (solid lines) and two with either a [100] screw dislocation core (dotted lines) or a [001] screw dislocation core (dashed lines) that vacancies can partition to. The systems with dislocation cores at some concentration develop excess Si vacancies-these Si vacancies are the one that partition to the dislocation core.                Table S1 Potentials used in our forcefields calculations. O (the normal oxygen in the crystal lattice) and OH (the oxygen in a hydroxyl group) have some unique but also some shared forcefields-O* represents both O and OH. The Morse potential for OH-H interactions was set to operate between 0 and 1.5 Å whereas the Buckingham potential for O-H interactions was set to operate between 1.5-10 Å. This ensures (with a sensible starting geometry) that the O-H bond is modelled by a Morse potential but the interaction of the oxygen in the OH group with the other Hydrogen in the vacancy is modelled with a Buckingham potential. If both potentials are set to operate from 0-10 Å then the hydrogen atoms either fall outside of the vacancy or into the centre of the vacancy (depending upon starting geometry) which does not match the more accurate predictions of DFT. [100](010)