Anisotropic thermal conductivity of antigorite along slab subduction impacts seismicity of intermediate-depth earthquakes

Double seismic zones (DSZs) are a feature of some subducting slabs, where intermediate-depth earthquakes (~70–300 km) align along two separate planes. The upper seismic plane is generally attributed to dehydration embrittlement, whereas mechanisms forming the lower seismic plane are still debated. Thermal conductivity of slab minerals is expected to control the temperature evolution of subducting slabs, and therefore their seismicity. However, effects of the potential anisotropic thermal conductivity of layered serpentine minerals with crystal preferred orientation on slab’s thermal evolution remain poorly understood. Here we measure the lattice thermal conductivity of antigorite, a hydrous serpentine mineral, along its crystallographic b- and c-axis at relevant high pressure-temperature conditions of subduction. We find that antigorite’s thermal conductivity along the c-axis is ~3–4 folds smaller than the b-axis. Our numerical models further reveal that when the low-thermal-conductivity c-axis is aligned normal to the slab dip, antigorite’s strongly anisotropic thermal conductivity enables heating at the top portion of the slab, facilitating dehydration embrittlement that causes the seismicity in the upper plane of DSZs. Potentially, the antigorite’s thermal insulating effect also hinders the dissipation of frictional heat inside shear zones, promoting thermal runaway along serpentinized faults that could trigger intermediate-depth earthquakes.

We designed a 2D numerical model to quantify the effect of antigorite's thermal conductivity anisotropy on the slab's thermal evolution.To perform the calculation, we employed a self-written MATLAB script.The slab was represented as a rectangular object with a thickness   , and a length   = 5  (Fig. S5).The slab thickness   was computed with the analytical solution reported by Ref 1 : where  is the slab age, here assumed as  ≈ 2.52 × 10 15  (80 ), and  is the thermal diffusivity set to  = 1 × 10 −6  2  −1 .The resulting slab was 120  thick, and 600  long.The value of 80  was chosen because it is considered as the upper limit for the slab thickening, representing a good proxy for a cold old slab 2 .
The cold slab subducted vertically into the mantle, i.e.,   = 90°, and it was progressively heated by the hot mantle (see Fig. S6).The slab was composed by three lithologies layered as follows: () a meta-basaltic crust on top (7  thick) 3 , () a layer of hydrated harzburgitic lithosphere in the middle (3  thick) 4 , and () a dry harzburgitic lithosphere keel at the bottom (110  thick) 5 .
The initial slab temperature profile at the trench was computed using the half-space cooling equation 1 : In this equation, we assumed a surface temperature of  0 = 300 , and an ambient mantle temperature of   = 1600 .The parameters  and  are the same of eq. 1, whereas  () represents the spatial coordinate.This equation produced a 1D temperature profile along the slab thickness, and the solution was truncated at 120  of depth.The temperature transect at the trench was then extended to the rest of the slab (Fig. S5).
Slab subduction was simulated from 120  to 230  depth by imposing a constant sinking velocity of   = 5   −1 .The entry depth of the slab in the mantle corresponds to the depth of the bottom part of the lithosphere, and it was chosen in order to avoid to simulate complex interactions between the slab and the overriding plate at the trench.The depth of The dominant heat transport mechanism inside the subducting slab is diffusion.The heat diffusion is described by the equation: where (/) indicates the temperature evolution over time (   −1 ), ∇ is the temperature gradient (   −1 ) in 3D space (  ,  ,  ), ∇ • is the divergence in 3D space (ʄ  /) + (ʄ  /) + (ʄ  /) ,  (  −3 ) is the density,  (   −1  −1 ) is the specific heat capacity and  (  −1  −1 ) is the thermal conductivity.In 2D the heat diffusion equation becomes: Given the model geometry of the sinking slab, the -direction represents the horizontal axis along the slab thickness, whereas the -direction represents the vertical axis along the slab length (Fig. S7).To investigate the anisotropy of antigorite's thermal conductivity we distinguished the horizontal   and vertical   components along the corresponding axis.
Antigorite was the only mineral in the model with   ≠   , whereas olivine   =   .
During the slab descent, the pressure and the temperature conditions were increased with depth.Pressure was increased by using the Preliminary Reference Earth Model (PREM) from Ref 7 (Fig. S8).The profile was built using a spline function with a 2 nd order polynomial and 10 intervals (Table S1 for the coefficients).

𝑇(𝐷
The compression-induced temperature increase ∆ was added to slab temperatures as well as the ambient mantle.The progressive increase of the  - conditions experienced by the descending of the slab determine the variation of the materials' parameters , ,  (see Note S3-S5).The governing equations of , , and  are reported in Note S3-S5 and Table S2-S3.

Note S2. Numerical Model
To solve the partial differential equation for 2D heat diffusion (eq.4.) we employed the finite difference (FD) method 9 .The domain was discretized into square cells with a constant grid space of ∆ = ∆ = 500 .The resulting grid mesh was composed by 242 × 1202 nodes.For the calculations we used the central difference method 9 , and the discretized heat diffusion equation reads: In this discretization, the temperature variations in a central node   is computed considering the four neighbouring nodes:  −1, (upper),  +1, (lower),  ,−1 (left), and  ,+1 (right).The To form the boundaries of the model, we prescribed four extra stencils of nodes in the four sides of the domain.Each node of the boundary was set to have isothermal conditions for the entire duration of the simulation.At depth above 120  (i.e., slab entry point in the mantle) the boundaries represented a 500--thick layer of ambient mantle.Therefore, the temperature of a given boundary node depended on the relative position of the slab along the adiabatic profile of the mantle (eq.6).Below 120  depth, the boundaries were designed to maintain the slab at the initial temperature conditions of the trench.To limit the internal heat diffusion, we set the left and right boundaries of the domain with the same temperature of the edges of the slab.Therefore, the right side of the domain, which represented the top of the vertically sinking slab, was set to 300 .The left side of domain representing the bottom of the slab was set to the same temperature of the truncated half-space cooling profile.The heat diffusion equation was solved implicitly with the direct solver algorithm 'mldivide', implemented in Matlab as the operator '\'.

Note S3. Thermal Conductivity 𝜦
Thermal conductivities of slab minerals are based on TDTR measurements at high pressures and variable temperatures.The thermal conductivities of different mantle phases, as a function of pressure   (), are expressed with a 4 ℎ order polynomial: The subscript '' indicates the given mineral phase (e.g.,   for olivine), the letters from - represents the various coefficients, and  0  stands for the thermal conductivity at ambient pressure.See Table S2 for a summary of all coefficients.
Temperature dependent   follows the relation reported in the present study and Ref 10,11 : where  298  is the thermal conductivity at high- and room- (298 ).The exponent   represents the -dependent coefficient for the given phase.
The aggregate thermal conductivity of each slab lithology was computed as the geometric average among the contributions of its constituent minerals 12 : where   is the thermal conductivity of the phase , and   is its volume fraction in the rock.
The lithospheric slab was assumed to be made entirely of olivine (  = 1), using the lattice thermal conductivity   reported by Ref 13 : Such assumption is justified by the argument that thermal conductivities of pyroxene and garnet, the main components of the meta-basaltic crust, are expected to be similar to the olivine 14,15 .In this model, we assumed a completely dry meta-basaltic crust and we concentrated all the hydrous minerals in 3--thick layer at the interface between the crust and the lithosphere (Fig. S6), as proposed by Ref 16 .The aggregate thermal conductivity of the hydrous layer   was computed by including the contributions of olivine  𝑂𝑙 13 and antigorite   measured in this study: A 15 % antigorite corresponds to an hydrous layer with 2% H 2 O 17,18 .To simulate different scenarios of lithosphere hydration, we changed the volume fraction of antigorite   in each model:   = 0.1 , 0.3 , 0.5 , and 1 .The 0.3 and 0.5 cases represent pervasive alteration scenarios, which can be found, for example, in the Atlantic oceanic lithosphere 19 .
The case of   = 1, on the other hand, represents an end-member scenario, which assumes the complete serpentinization of the hydrous layer.
For our present study, we run 10 different models (Figure S9a-c indicates the fundamental cell units of a given mineral, and 6.022 is linked to Avogadro number (   = 6.022 × 10 23 ).To compute  -dependent density we took the thermal expansion coefficient  of each mineral from literature.  (, ) =  0  +    +    +    + ℎ   2 (olivine) The reduction of cell volume caused by temperature increase ∆ is described by the relation: Density was finally obtained by using the molar weight    (  −1 ) of each mineral: Table S3 summarizes all coefficients.Data for olivine were taken from:    () (Downs et al., 1996),   (, ) 20 .Antigorite data were taken from:    () 21 ,  𝐴𝑡𝑔 22 .
We computed the density of each lithology as the weighted average of the constituting minerals (Figure S9d): where   indicates the volume fraction of antigorite in the rock.

Note S5. Specific Heat Capacity 𝑪𝒑
The specific heat capacity   (   −1  −1 ) of each mineral phase was obtained from the measurements of volumetric    (   −3  −1 ) and molar    (   −1  −1 ) heat capacities, see Table S3 for a summary of all coefficients.We extrapolated olivine's specific heat capacity   from the dataset of Ref 20 : We extrapolated the specific heat capacity of antigorite   from the dataset of Ref 22 :   (, ) =    6 +    5 +    4 +    3 +    2 +    1 −  0  (31) We computed the specific heat capacity of each lithology as the weighted average of the constituting minerals (Figure S9e): (33)  where   indicates the volume fraction of antigorite in the rock.

Note S6. Code limitations:
Our model setting presents several limitations due the assumptions employed in the calculations: (1) We computed the heat transfer equation for 2D case.The resulting slab can only be heated from four sides (left, right, bottom, top), and thus is expected to be colder than the 3D geometry (width is missing).
(2) We used a simplified petrology for the three lithologies of the slab, by assuming to be exclusively composed by olivine and antigorite.The aggregate thermal conductivity of each lithology should be calculated including the contribution of other major phases: clinopyroxenes and orthopyroxenes, garnets (basaltic crust and oceanic lithosphere); Al-rich phases (phlogopite, phengite, chlorite) and SiO2-polymorphs (sediment blanket), other hydrous phases (e.g., talc, brucite, humite group).However, it should be considered that the lack of thermal conductivity datasets on major mantle minerals has been progressively filled by our recent measurements.It's only a matter of time before more comprehensive aggregate  estimates of mantle rocks will be available.
(3) The phase proportions in each lithology should vary accordingly to their - stability field.In our calculation we simply assumed that antigorite remains stable at any P-T condition.Future models should compute the breakdown of antigorite as the temperature increase.
(4) Our model is purely thermal, and subduction occurs with a constant sinking velocity.
Future models should be thermo-mechanical, thus including buoyancy-driven subduction, and the variable sinking velocity should arise self-consistently from the interactions between the slab and the mantle.
(5) The contribution of additional heating source should be considered, e.g.radioactive elements (enriched in the basaltic crust and in the sediments), frictional heating, and positive/negative latent heat at the phase transitions.(6) The domain in our model is limited to the slab, and a 500 m thick surrounding mantle.
In nature, the slab heating causes the corresponding cooling of the mantle around the slab.Therefore, the temperature of the ambient mantle changes as the slab heats up.
Future models should include larger volumes of the ambient mantle to reproduce this effect.(7) In our model subduction starts at 120 km to avoid the interactions between the slab and the overriding plate.These interactions, however, influence the temperature field of the slab in the shallow portion of the subduction zone, thus influencing its P-T-t path before being fully decoupled from the overriding plate.Future models should include these interactions.
(8) Corner-flow models 16 show that as the slab subducts there is an upwelling of hot material from deep mantle regions toward the mantle wedge.Moreover, part of the shallow cold mantle is dragged at depth by the subducting slab, thus creating a cold nose in the mantle wedge.This cold material can protect the slab from the high temperatures of the deep mantle, thus maintaining the slab surface cold for a longer time.Future models should also include the corner flow in the mantle wedge.

Note S7. Code benchmark:
We verified the reliability of the numerical model by comparing it against a known benchmark case.We designed a 2D model with Gaussian temperature distribution.Benchmark settings are summarized in Table S7.We solved the Fourier's heat diffusion law numerically and benchmark it against the analytical solution for the same problem setting (Fig. S10): where  is the temperature as a function of position  and time ,  max is the peak temperature in the Gaussian profile,  is the thermal diffusivity (/( * )), and  the amplitude of the Gaussian distribution.To quantify the goodness of our discretization we computed the misfit between analytical and numerical solution, at a given node position.We further computed the  between the analytical   and numerical solution   , and the  2 norm (Fig. S11): where  is the total length of the domain and ∆ the grid spacing.As reported in Fig. S10, the  2 << 1 , thus indicating a negligible misfit.The benchmarking proves the chosen discretization is sufficient to solve numerically the heat diffusion equation.

Note S8. Assumption of perfect alignment of single-crystalline antigorite along the slab dip angle:
Several field observations, deformation experiments, and theoretical calculations have reported that antigorite has a strong crystal preferred orientation (CPO), in which its [001]   direction tends to orient perpendicularly to the shear direction.For example, Padrón-Navarta et al. 23 reported a strong foliation in natural serpentinites samples (Cerro del Almirez, Spain): antigorite's [001] shows a peak of relative density of 4.9-17.4m.u.d.(multiples of uniform distribution) along the Z-axis, which is perpendicular to the shear direction (X-axis).Moreover, the majority of antigorite grains orient their [001] direction ±30º from the Z-axis.Similar behaviors have also been reported by Nishii et al. 24 (4.7-18 m.u.d.) and Jung 25 (5 m.u.d.).On these bases, it is reasonable to assume that, during subduction, antigorite would orient its [001] direction normal to the main shear direction, i.e., along the slab dip.
Unfortunately, it is not possible to find ophiolites of subducted slab which maintained the CPO gained during subduction.Moreover, to our knowledge, it is not yet possible to infer the orientation of antigorite inside the slab from seismology.Satta et al. 26 reported that in order to see 1 second delay on the shear-velocity, the seismic wave has to cross a >10 km-long path of fully serpentinized rocks (100% alteration), and the incident angle of the seismic wave must be parallel to the foliation plane.Given such conditions and the limited time sensitivity of modern seismometers (~ 10 milliseconds), it is almost impossible to tell the orientation of a partially serpentinized rock layer (<70% alteration) (see Allen et al. 27 and Cooper et al. 28 ) which is crossed by a non-ideally oriented seismic wave for less than 10 km.Therefore, it is reasonable to assume that antigorite's [001] direction is oriented perpendicularly to the slab dip within ±30º.
Moreover, our TDTR measurements report the upper and lower bounds of antigorite's thermal conductivity depending on the orientation of the crystal.We can define  as the incident angle between the heat flux and the foliation plane of antigorite.On one hand, if the [001] direction is parallel to the heat flux (i.e., perpendicular to the {001} foliation planes  = 90º), the heat will flow slowly because of the low Λ 001 .On the other hand, if the [001] direction is perpendicular to the heat flux (i.e., parallel to the {001} foliation planes  = 0º), the heat will flow faster because of the high Λ 010 .If the antigorite is oriented with an incident angle between 0º <  < 180º and  ≠90º, we can assume that the thermal conductivity would be a value between Λ 001 and Λ 010 .In this case, we can use a simple function: Here we used a geometrical average because Λ 001 is much lower than Λ 010 , and hence it will act as a bottleneck to the lattice heat transport.With this equation we can calculate antigorite's thermal conductivity as a function of its orientation with respect to the heat flux (Table S8, Fig. S13).
Equation 36 is a fair approximation of the thermal conductivity of antigorite when 0º ≤  ≤ 180º:    remains low for high incidence angle (45º ≤  ≤ 135º) because heat would flow through several foliation planes.However,    is high when the incidence angle is low (  ≤ 30º or  ≥ 150º) because heat would flow on the foliation planes, see Fig. S13.
Heat flux is always oriented along the direction of the maximum temperature gradient, which, in our case, is perpendicular to the slab surface (Fig. S6).Therefore, we can assume that the incident angle between heat flux and antigorite's foliation planes is high:  = 90º±30º.This range of orientations does not show a major difference in antigorite's thermal conductivity: ~20% at ambient pressure which becomes less than 10% at high pressure (see the plateau of  001 in Fig. 1 of the main text).Even at  = 30º and  = 120º, the  010 / 001 anisotropy remain large: ~3.4 at ambient pressure, and ~2.3 at P=7 GPa.For these reasons, we believe that the assumptions made to design our simplified model of heat flow through an anisotropic medium are sufficient to describe the physics of the phenomenon, and our conclusions are still valid even if antigorite is not perfectly aligned.

Note S9. Potential effects of grain boundary:
As we show in the main text, our single-crystal data along the b-and c-axis provide important basis that brackets the thermal conductivity of polycrystals.Its applicability to estimate the "average" thermal conductivity of a rock that likely contains polycrystals, depends on many factors, such as the realistic fraction of crystals along b-and c-axis.This, in turn, is influenced by: 1) the grain size; 2) the misorientation of the polycrystalline aggregate with respect to the local stress state; 3) the presence of other phases in the aggregate; 4) the largescale stress field imposed by the tectonic settings.
As we discuss in Note S8, it has been reported that the maximum orientation density of the [001] direction along the Z-axis is ~4.7-18 m.u.d., and most of the grains orient their [001] within a ±30º angle with respect to the Z-axis.The thermal conductivity difference between a fully aligned aggregate and a slightly misaligned one is ~8-20%, whereas its anisotropy remains high (~2.3-3.4).In order to compute the aggregate thermal conductivity of an antigorite mineral assemblage, one has to determine: 1) the distribution density of the [001]   orientation with respect to the Z-axis (which should be less than the ~4.7-18 m.u.d peak); and 2) the grain size of each crystal.With this information it is then possible to compute a weighted average considering what percentage of the bulk volume is oriented in which direction.There are codes (e.g., D-Rex 29 ), which allow to compute the orientation and the grain size of crystals produced by an external stress field.These codes, however, are limited to olivine and enstatite slip systems.Implementing antigorite slip systems in these codes to produce CPO pole figures is beyond the scope of this paper.Our simplified approach is sufficient to prove the importance of antigorite anisotropy on slab's thermal evolution.
The presence of misaligned grain boundaries can generate additional scattering of the heat-carrying phonons and result in a thermal resistance, i.e., thermal conductivity reduction 30 .
The scattering due to tilt boundaries, however, seems to be minor with respect to scattering in the inter-grain regions, which is where the crystallographic disorder concentrates 30 .Moreover, it has been shown that the thermal resistance across the basal plane of mica (a phyllosilicate like antigorite) decreases as the planes are forced together by the increasing pressure 31 .As suggested by Wood 32 , the fall in mica thermal conductivity reported by Powell and Griffiths 31 is attributed to the grain size reduction due to the crushing of the sample at high pressure.From these experimental evidences and theoretical considerations, it seems that phonon scattering at the grain boundaries of a phyllosilicates aggregate is mostly affected by the grain size rather than the misalignment of the crystals.
We therefore expect that the misalignment of antigorite crystals in a serpentinized rock has a minor impact on the aggregate thermal conductivity of the mineral assemblage, because 1) the misalignment is within a narrow angle (±30º) and 2) the misalignment tends to reduce at high pressure and under high shear stress.
It will be interesting to investigate the effects of grain size on antigorite thermal conductivity.Padrón-Navarta et al. 23 report a natural sample of serpentinite with a grain size of few tens of μm, which is at the same order of magnitude of the single crystal samples we measured.If the reduction of thermal conductivity due to smaller grain size is proved experimentally in the future, it would further strengthen our point: a fine-grain aggregate of aligned antigorite crystals produced by an shear stress field would create an even better thermally insulating layer.Thus, the antigorite's anisotropic thermal conductivity is expected to play an important role in influencing the temperature evolution of the slab down to ~250 km depth.S6, ΛAtg=3.7 W m -1 K -1 (red curve) offers a best-fit to the raw spectrum.The ratio -Vin /Vout is most sensitive to antigorite's thermal conductivity at few hundred picosecond (ps) delay time 33,34 .A 10% test error for ΛAtg (green and blue curves) causes the fitting curve deviating from the data.Such high sensitivity demonstrates that based on the high-quality data, our thermal model fitting and the obtained ΛAtg are precise and reliable.We simplified the slab as a 120 × 600  rectangle, and we discretized the domain using square grid cells.Each cell is bounded by 4 nodes at the corners with a constant grid spacing of 500 .We used the [, ] coordinate system to identify each node:  identifies the rows (-direction) and  identifies the columns (-direction).We used a Scientific Colour Map to plot the temperature inside the slab (see caption of Fig. S5).    >  10 001 (Fig. 5a) at shallow depth due to the low thermal conductivity, but it becomes  10 010 <  10 001 at greater depth because the hydrous layer is more thermally insulated than the   010 configuration.Compared to antigorite, the heat flux through the hydrous layer in the dry olivine reference ( 7  , solid olive-green line in Fig. 5b) is higher than  7 001 (horizontal insulator), but lower than  7 010 (vertical insulator), because the thermal conductivity of antigorite along both direction is  010 >   >  001 .S8).We used our TDTR measurements to set the two endmembers values of antigorite's thermal conductivity: cross-plane  001 (min); and in-plane  010 (max).The trend reproduced by eq.36 (red line) can be physically justified by the following geometric considerations: 1) Λ  is minimum at  = 90°, because the heat flows through all foliation planes, i.e. parallel to the [001] direction    =  001 .Note that, for high incidence angles ( 60°<  < 120°), a small misalignment of  = ±30° results in    ~001 + 0.25   −1  −1 (red area). between the heat flow and the foliation planes {001}.Each thermal conductivity was computed using eq.36:    () = (Λ 001 )  * (Λ 010 ) 1− .This empirical equation can be use to compute    of non ideally oriented antigorites crystals.

Figure S1 .
Figure S1.Representative (a) X-ray diffraction and (b) Raman spectrum of our natural antigorite sample.In (a), the blue curve is our measurement result and the red peaks are from antigorite's X-ray diffraction PDF database.In (b), the Raman intensity is in arbitrary unit (a.u.).

Figure S2 .
Figure S2.Schematic drawing of the high-pressure TDTR experimental setup in a DAC at room temperature.A single-crystal antigorite coated with Al film was placed in the DAC chamber.A ruby ball and silicone oil were used as the pressure calibrant and the pressure transmitting medium, respectively.The thermal conductivity of antigorite was measured by optical pump-probe method (see the main text for the details).

Figure S3 .
Figure S3.Example TDTR spectrum (open circles) for antigorite at 8 GPa and room temperature along c-axis loaded with Ar as the pressure medium.The data are fitted by thermal model calculations (color solid curves).Using the input parameters listed in TableS6, ΛAtg=3.7 W m -1 K -1 (red curve) offers a best-fit to the raw spectrum.The ratio -Vin /Vout is most sensitive to antigorite's thermal conductivity at few hundred picosecond (ps) delay time33,34 .A 10% test error for ΛAtg (green and blue curves) causes the fitting curve deviating from the data.Such high sensitivity demonstrates that based on the high-quality data, our thermal model fitting and the obtained ΛAtg are precise and reliable.

Figure S4 .
Figure S4.Sensitivity tests of the thermal model calculations to input parameters for antigorite at 8 GPa and room temperature loaded with Ar as the pressure medium.Here antigorite's thermal conductivity, ΛAtg, is fixed at 3.7 W m -1 K -1 , as derived from Fig. S3, using parameters listed in Table S6.(a) and (b) Even if the thicknesses of Ar (hAr) and antigorite (hAtg) change

Figure S5 .
Figure S5.Schematics of the slab used in our models.The slab was simplified as a 120 × 600  rectangle.The colour map represents the initial temperature of the slab: dim colours represent areas with low  , whereas bright colours indicate high  .To plot temperature, we used the Scientific Colour Maps 'lajolla'.

Figure S6 .
Figure S6.Schematics of slab lithologies.We subdivided the slab into different domains: oceanic crust (7--thick), hydrous layer (3--thick), and lithospheric slab (110--thick).Each lithology is characterized by a different mineral assemblage (see main text).We added an extra stencil of nodes at the four sides of the slab to represent the ambient environment (surface lithosphere or deep mantle).The red arrows indicate the heat flow from the hot ambient mantle toward the cold slab.

Figure S7 .
Figure S7.Schematics of discretization of the slab.We simplified the slab as a 120 × 600  rectangle, and we discretized the domain using square grid cells.Each cell is bounded by 4 nodes at the corners with a constant grid spacing of 500 .We used the [, ] coordinate system to identify each node:  identifies the rows (-direction) and  identifies the columns (-direction).We used a Scientific Colour Map to plot the temperature inside the slab (see caption of Fig.S5).

Figure S8 .
Figure S8.(a) Pressure profile of the mantle 7 .(b) Adiabatic temperature profile in the mantle 8 .(c) Temperature gradient in the mantle 8 .

Figure S10 .
Figure S10.Benchmark: comparison between analytical and numerical solutions.The two subplots illustrate two perpendicular transect of the 2D Gaussian distribution of temperature: (a) along x-direction; (b) along y-direction.In the figure are shown 5 different snapshots of the temperature evolution (0;0.25;0.5;0.5;1Myrs).We used a Scientific Colour Map (see caption ofFig.S5) to plot the progressive cooling of the domain from bright (hot) to dim (cold).The solid lines represent of the analytical solutions, whereas the e markers indicate the value of the numerical solution at the corresponding node (every 50th nodal point is plotted here).

Figure S11 .
Figure S11.Benchmark: (a)  2 norm evolution over 1 Myrs of simulated time, (b) misfit between the analytical and numerical solution.

Figure S12 .
Figure S12.Radial profile of heat flux  (  − ) computed at both sides of hydrous layer.(a) Internal side (10  from slab surface).(b) external side (7  from slab surface).In each model, we used the Fourier Law of heat conduction to compute the heat flux:  = −∇.Each line indicates the  profile of one model: green solid lines for set 1 as references, dotted purple lines for set 2 with vertical insulator, and dashed teal lines for set 3 with horizontal insulator.The two boxes below represent the two configurations of antigorite used in the model (see Fig.3 in the main manuscript), while the vertical black dashed lines indicate the two sides of the hydrous layer.The heat flux entering the hydrous layer from the external side ( 7 , Fig. 5b) is higher in the model set 2 (vertical insulator, purple dotted lines) compared to model set 3 (horizontal insulator, teal dashed lines),  7 010 >  7 001 .In the   001 configuration,  10010 >  10 001 (Fig.5a) at shallow depth due to the low thermal conductivity, but it becomes  10 010 <  10 001 at greater depth because the hydrous layer is more thermally insulated than the   010 configuration.Compared to antigorite, the heat flux through the hydrous layer in the

Table S1 .
2)    is maximum at  = 0° and  = 180°, because the heat flows on the foliation plane, i.e. parallel to the [010] direction    =  010 .In this case, at low incidence angles ( 0°<  < 30° and 150°<  < 180°), a small misalignment of  = ±30° results in    ~010 − 2   −1  −1 .In this figure are reported only the values at ambient ,.The layered structure inside the plot represents the heat flowing through an antigorite crystal: the orange arrows indicates the heat flow with different incidence angles ( = 90°, 45°, 0°); the purple areas represent antigorite planes characterised by high thermal conductivity  010 ; the teal areas represent the cross-plane regions characterized by low thermal conductivity  001 .Coefficients to compute the mantle pressure profile.

Table S3 .
Coefficients to compute unit cell volume, thermal expansion coefficients, and heat capacity.

Table S7 .
Parameters used to compute the analytical solution for a Gaussian distribution of the temperature field

Table S8 .
Antigorite thermal conductivity    , computed for different incidence angles