高压下Py-FeO2、Py-FeOOH 和ε-FeOOH 晶体结构与弹性特征的第一性原理研究

顾小雨，刘 雷

（中国地震局地震预测研究所高压物理与地震科技联合实验室，北京 100089）

地球表面70%被水覆盖，同时更大量的水以不同形式赋存在地球内部。地表的水可以通过板块俯冲作用被带到地球内部，随着深度的增加，地球内部的水将经历复杂的迁移过程[1]。水在地球深部中的运移和储存会对地球内部动力学和物质循环演化产生影响，同时还影响地壳和水圈的演化规律。水的存在显著影响成岩矿物的化学和物理性质，包括矿物相关系、熔化温度、流变特性、电导率和地震速度等。在地幔深度范围内，水主要以羟基形式存在于矿物晶格[2]。最近对含水矿物超石英（stishovite）的原位激光与X 射线衍射（X-ray diffraction，XRD）实验表明，超石英可以在地幔温度压强条件下稳定存在，超石英有可能成为下地幔中运输和储存水的关键矿物[3]。Li 等[4]使用从头算分子动力学和热力学方法计算了20～135 GPa 和5 800～5 000 K 条件下铁与硅酸盐熔体之间的水的分配比，验证了水的亲铁性，发现地核可能是一个包含地球大部分水的大型水库，这可能是引起地心低密度的原因。但是，目前水如何输送到地球的深部，以及在哪些阶段进行，仍然存在很大的争议。地幔物质的性质或许可以通过某些地震学观测约束得到，但是理解地幔的具体组成需要完备的、高质量的高温高压矿物数据，如密度、稳定性、弹性性质[2-27]等。因此，研究高温高压下含水矿物的性质，对了解地球深部的结构、矿物组成以及深部水循环具有重要意义。

本研究选择地球深部重要的含水相矿物黄铁矿型FeO2（Py-FeO2）、FeOOH（Py-FeOOH）以及ε-FeOOH[5-7]作为研究对象。Py-FeO2和Py-FeOOH 是2016 年在高温高压实验中被发现稳定存在的新矿物[5-6]，实验和理论结果都证明了从ε-FeOOH 到黄铁矿型FeOOH 的转化关系，ε-FeOOH 和黄铁矿型FeOOH 之间的相边界位于60 GPa、1 500 K[7]。

Hu 等[5]使用第一性原理计算和实验，首次确定了76 GPa、1 800 K 条件下稳定存在的P 相（黄铁矿型FeO2）。Bo 等[8]利用密度泛函理论（density functional theory，DFT）和DFT + DMFT 方法研究了FeO2的电子结构和结构性质，FeO2的晶体结构由Fe2+和O22-二聚体组成，其中Fe 离子被八面体O 原子包围，FeO2在高压下表现出金属-绝缘体（MIT）转变（由氧二聚体键长度变化驱动的绝缘类型）。Hu 等[6]进一步确定了P 相至少稳定到133 GPa、2 000 K 条件下，实验中还发现了具有黄铁矿型结构的FeOOHx，利用第一性原理计算方法约束和量化了P 相（Py-FeOOHx）中的氢含量（0.39 ＜x＜ 0.81）。Nishi 等[7]发现了新黄铁矿型FeOOH 相，可在地幔底部的条件下稳定存在，并最终分离成Fe2O3和H2O。Zhang 等[9]进行了第一性原理计算，评估了黄铁矿型FeO2的稳定性和弹性性质。Huang 等[10]进行了第一性原理计算，发现氢的加入显著影响了高压下黄铁矿型FeO2H-FeO2体系的结构特性、压缩行为和各向异性。Liu 等[11]利用高压实验和理论计算分析了Py-FeOOHx的声速结果，发现铁和水能在86 GPa、2 200 K 条件下反应生成Py-FeOOHx。Yuan 等[12]基于原位XRD 实验和第一性原理计算，研究了Fe-H2O 在高压下的相变特征，发现在低于78 GPa 的条件下产生FeHx和FeO，高于78 GPa 的条件下生成Py-FeOOHx和FeHx。Zhuang 等[13]利用Kawai 型多砧压机和XRD 实验证明了在Py-FeOOHx转化为ε-FeOOH 的过程中，ε-FeOOH 发生了部分脱氢反应，以ε-FeOOHx的形式存在于下地幔中。Deng 等[14]利用第一性原理分子动力学研究了P 相的固液两相的热弹性性质，发现通过相对少量的熔融FeOOHx可以解释观察到的地幔边界上方超低速区（ultra-low velocity zones，ULVZs）的地震异常。Hou 等[15]结合高压电导率实验、拉曼光谱测量和第一性原理模拟，研究了黄铁矿型FeOOHx（x≤ 1）中氢的状态，发现当压强超过73 GPa时，对称的羟基键被软化，H（或质子）在其晶体位置附近变得分散。高度扩散的氢为电荷和质量提供了新的输运机制，这决定了地球深部的电导率和磁性等地球物理性质，以及地幔底部中氧化还原、氢循环和氢同位素混合的地球化学过程。Thompson 等[16]利用第一性原理计算了黄铁矿体系FeOOHAlOOH 在60～140 GPa 条件下的弹性性质和地震波速，发现黄铁矿型(Al0.5～0.6, Fe0.4～0.5)OOH 可以很好地吻合下地幔中超低速区的压缩和剪切速度，富含铝的CaCl2型(Al,Fe)OOH 可能形成大的低剪切速度区域（large low shear-wave velocity provinces，LLSVPs）。

Bendeliani 等[17]在9.0 GPa、400 K 条件下合成了ε-FeOOH，XRD 和中子衍射实验[18-20]表明，ε-FeOOH 具有空间群Pmn21的正交结构。实验发现[21-22]，α-FeOOH 在5 GPa、200 ℃的条件下需要较长的反应时间才能生成ε-FeOOH。Suzuki[23-24]利用同步加速XRD 实验，得到了8.6 和11 GPa、700 K 温压条件下ε-FeOOH 的晶体结构参数和体积模量（K）。Otte 等[25]利用DFT 确定了ε-FeOOH 的稳定性，同时得到了ε-FeOOH 的结构、磁性和电子性。Gleason 等[26]采用实验与DFT 计算相结合的方法，研究了ε-FeOOH 的结构和电子自旋跃迁。Thompson 等[27-28]使用DFT 研究了ε-FeOOH 的结构和弹性性质，并利用XRD 实验对ε-FeOOH 的高压行为进行研究，发现在0～75 GPa 压力范围内ε-FeOOH 有2 个连续的结构相变，(18 ± 1) GPa 时由于压力诱导的氢键对称化而产生的P21nm→Pnnm结构变化，在(45 ± 2) GPa时经历了等结构的一级相变。

Py-FeO2、Py-FeOOH 和ε-FeOOH 3 种矿物作为Fe-H2O 体系尚未被系统地研究，Py-FeOOH 和Py-FeO2的弹性性质还有待进一步探索。第一性原理方法在探索物质在高温高压下的结构、弹性、谱学等方面取得了丰硕的成果，且计算精度不断提高[29-39]。为此，本研究通过第一性原理计算分析0～350 GPa 范围内的Py-FeO2和Py-FeOOH 以及0～170 GPa 范围内的ε-FeOOH 的晶体结构、弹性模量以及波速等性质，系统地比较3 种矿物的物性差异，为了解高压下Fe-H2O 体系矿物物性演化及其在地球深部赋存特征提供数据支持。

1 方法与理论

采用基于DFT 的第一性原理方法，使用Materials Studio-CASTEP[40]模块，计算了Py-FeO2、Py-FeOOH 和ε-FeOOH 在不同压强下的结构变化。晶胞结构优化计算使用基于PBE（Perdew-Burke-Ernzerh）交换关联函数的广义梯度近似方法（generalized gradient approximation，GGA），GGA 方法更适用于模拟含氢键体系矿物[24,27,35]。在收敛性测试的基础上，平面波基组展开的截断能均为571.4 eV，布里渊区设置K点的网格参数为6 × 6 × 6（Py-FeO2）、6 × 6 × 6（Py-FeOOH）、5 × 6 × 8（ε-FeOOH）。

第一性原理计算的第一步需要对Py-FeO2、Py-FeOOH 和ε-FeOOH 的晶体模型进行几何优化，求得局域最稳定结构，晶体结构模型如图1 所示。Py-FeO2和Py-FeOOH 属于立方晶系（空间群Pa3，原子序数Z= 4）[5-7]，ε-FeOOH 为正交晶系（空间群Pmn21，Z= 2），结构中的Fe 和O 原子形成了CaCl2型结构[25]。不同压强下ε-FeOOH、Py-FeOOH 和Py-FeO2的晶胞参数及结构内各原子位置都通过BFGS 方法进行了优化。通过应力-应变关系计算弹性常数，控制应变大小为0.05，以保障ε-FeOOH、Py-FeOOH 和Py-FeO2处于弹性形变状态。弹性模量利用Voigt-Reuss-Hill 近似方法计算得出。

图1 Py-FeO2、Py-FeOOH 和ε-FeOOH 的晶体结构Fig. 1 Crystal structure of Py-FeO2, Py-FeOOH and ε-FeOOH

2 结果与讨论

首先对Py-FeO2、Py-FeOOH 和ε-FeOOH 的晶体结构进行优化，得到稳定结构的晶格参数，见表1。Py-FeO2的晶格常数和晶胞体积与前人计算值[5]的相对偏差分别为0.14%和0.40%，与实验值[5]的相对偏差分别为2.4%和1.7%。Py-FeOOH 的晶胞体积与前人计算值[7]的相对偏差为2.3%，晶格常数和晶胞体积与实验值[5,7]的相对偏差都不超过0.5%。ε-FeOOH 的晶胞体积与前人计算值[25]的相对偏差为0.02%，晶胞体积与实验值[17,23-24]的最大相对偏差为14%，晶格常数a、b、c与实验值[17,23-24]的最大相对偏差分别为4.8%、4.7%、3.4%，前人研究结果中，ε-FeOOH 的体积模量计算值[24]与实验值[17,23-24]的最大相对偏差为32%。ε-FeOOH 的晶格常数、晶胞体积和体积模量与实验值的偏差略大，主要是由于计算中使用了GGA 方法，从而低估了结构内原子间的内聚力[17,23-24,41-43]。总体而言，本研究的计算结果与文献中的计算和实验结果[5-7,17,23-25]的一致性较好，表明采用的计算方法和参数设置是合理可靠的。

表1 Py-FeO2、Py-FeOOH 和ε-FeOOH 的晶格参数和体积模量Table 1 Lattice parameters and bulk moduli of Py-FeO2, Py-FeOOH and ε-FeOOH

2.1 晶体结构和密度变化

0～350 GPa 条件下的Py-FeO2和Py-FeOOH 以及0～170 GPa 条件下的ε-FeOOH 的晶格常数如图2 所示。Py-FeO2和Py-FeOOH 的结构为立方晶系（a=b=c），晶格常数随压强的增加呈逐渐减小的趋势，晶胞体积也随着压强的增加而降低。ε-FeOOH 的结构为正交晶系（a、b、c不相等），晶格常数随着压强的增大而逐渐减小，0～170 GPa区间内c轴长度相比a、b轴长度大幅减小，但在33 GPa 附近有异常突变，a、b轴长度突然增加，c轴长度突然减小。c轴上只有共边的氧原子和氢键，相比有共边的铁原子以及铁氧键的a、b轴，更容易被压缩。本研究在确定计算方法和参数后，具体计算了0～350 GPa 压力条件下Py-FeO2、Py-FeOOH 以及0～170 GPa 压力条件下ε-FeOOH的晶格常数和密度。

图2 Py-FeO2、Py-FeOOH 和ε-FeOOH 的晶格常数随压强的变化Fig. 2 Lattice parameters of Py-FeO2, Py-FeOOH and ε-FeOOH under high pressures

ε-FeOOH 的晶格常数在33 GPa 附近发生了突变，可能是ε-FeOOH 结构在压强驱动下出现了氢键对称化现象。氢键对称化与压强关系的相关研究表明，氢键对称化现象会在30～40 GPa 的压强驱动下在各种含有弱氢键的矿物相中产生，如冰的Ⅶ相、D 相、In-OOH 相以及含水硫化物系统等[44-47]。矿物中的羟基通常以不对称结构的O―H···O 开始，氢原子通过强短共价键（O―H）与一个氧原子结合，氢原子与另一侧的氧原子通过弱长氢键（H···O）结合。随着深度和压强的增加，羟基的共价键（O―H）延长，弱氢键（H···O）缩短，直到两个键的长度相等，H 到达两个相邻氧原子的对称中心。物质晶体结构在高压下发生氢键对称化现象，在具体数据中表现为羟基的共价键（O―H）和弱氢键（H···O）的键长变为O···O 键距离的一半[44]。

在33 GPa 以上，ε-FeOOH 结构中羟基的键长与O···O 键之间距离的一半几乎相等，如图3 所示。在ε-FeOOH 相这类含有弱氢键的体系中，常压下O···O 之间的距离一般为2.6 Å，羟基的键长一般为1.0 Å；随着压强的增加，ε-FeOOH 的晶胞体积逐渐缩小，33 GPa 时，O···O 之间的距离逐渐缩小至2.338 Å；羟基的键长随着压强的增加反而逐渐延长，33 GPa时羟基的键长逐渐增加至1.170 Å，此时羟基键长为O···O 距离的一半；随着压强继续增加（33 GPa之后），ε-FeOOH 结构中O···O 之间的距离和羟基键长继续减小，羟基键长与O···O 距离的一半相等。

图3 ε-FeOOH 结构中羟基键长和氧氧键长随压强的变化Fig. 3 Hydroxyl bond length and oxygen bond distance of ε-FeOOH under high pressures

P y-F e O O H 的结构中存在羟基，但P y-FeOOH 的晶格常数在30～40 GPa 没有发生突变，原因是Py-FeOOH 晶体结构内不存在不对称的O―H···O，而且在前人的实验[6-7]中Py-FeOOH 的压强（大于70 GPa）超过发生氢键对称化的压强条件（30～40 GPa）。Py-FeO2的结构内不存在羟基，结构中O―O 键长随压强变化的趋势也与ε-FeOOH 结构不同，计算发现结构中O―O 键长随压强的增加逐渐减小（见图4），与晶格常数随压强变化的趋势相同。

图4 Py-FeO2 结构中氧氧键长随压强的变化Fig. 4 Oxygen bond length of Py-FeO2 under high pressures

Py-FeO2、Py-FeOOH 和ε-FeOOH 的晶胞密度随压强的增加而增加（见图5），Py-FeO2、Py-FeOOH 和ε-FeOOH 的密度在同等压强条件下依次减小，Py-FeO2的密度最高，Py-FeOOH 的密度比Py-FeO2降低9.4%，ε-FeOOH 的密度比Py-FeO2降低14.4%。计算得到的密度随压强的变化趋势与实验结果[5,7]基本相符，Py-FeO2的晶胞密度与实验结果[5]的相对偏差为0.5%，Py-FeOOH的晶胞密度与实验结果[7]的相对偏差为0.8%，ε-FeOOH 的晶胞密度与实验结果[7]的相对偏差为2.5%。

图5 Py-FeO2、Py-FeOOH 和ε-FeOOH 的密度随深度的变化Fig. 5 Variation of density of Py-FeO2, Py-FeOOH and ε-FeOOH with depth

在0～3 000 km 深度范围内，Py-FeO2、Py-FeOOH 和ε-FeOOH 的密度高于典型地球密度结构模型（见图5），与地幔的平均密度（5 g/cm3）相比[11]，Py-FeO2、Py-FeOOH 和ε-FeOOH 的密度更高。Py-FeO2、Py-FeOOH 和ε-FeOOH 的密度曲线在2 898 km（核幔边界）附近与 PREM 和 AK135 模型相交，地球核部密度远高于Py-FeO2、Py-FeOOH 和ε-FeOOH。与周围的地幔相比，由于其明显更高的密度，Py-FeO2、Py-FeOOH 将继续向地球深部下沉。研究发现，在超过100 GPa 的压强条件下，Py-FeO2、Py-FeOOH 会与AlOOH 和MgSiO4H2形成固溶体[7]，Tsuchiya 等[48]预测AlOOH 能以黄铁矿型结构稳定存在至170 GPa，形成的固溶体会降低整个体系的密度，使之更接近PREM 模型。因此，Py-FeO2、Py-FeOOH 和ε-FeOOH 及固溶体的存在可能会对地幔的密度产生重要影响。

2.2 弹性性质

矿物的弹性性质是建立地球内部模型的重要参数，弹性性质包括弹性常数和弹性模量。为了研究Py-FeO2、Py-FeOOH 和ε-FeOOH 的弹性性质，分别计算了10～350 GPa 范围内Py-FeO2和Py-FeOOH 以及0～170 GPa 范围内ε-FeOOH 的弹性常数，结果见表1。

矿物的稳定性不仅包括能量的稳定性，还包括力学稳定性，力学稳定性可以采用Born-Huang 稳定性判据[49]进行判断。

(1) Py-FeO2和Py-FeOOH 的结构是立方晶系，有3 个独立弹性常数（C11、C12、C44），稳定性条件为：Cii＞ 0（i= 1, 4），C12＞ 0，C11-C12＞ 0，C11+ 2C12＞ 0。

(2) ε-FeOOH 的结构是正交晶系，有9 个独立弹性常数（C11、C12、C13、C22、C23、C33、C44、C55、C66），稳定性条件为：Cii＞ 0（i= 1, 2, 3, 4, 5, 6），C11+C22+C33+ 2（C12+C13+C23） ＞ 0，C11+C22- 2C12＞ 0，C11+C33-2C13＞ 0，C22+C33- 2C33＞ 0。

Py-FeO2和Py-FeOOH 在0～350 GPa 内满足立方晶系结构的稳定性条件，ε-FeOOH 在13～170 GPa范围内满足正交晶系结构的稳定性条件，表明它们均具有力学稳定结构。ε-FeOOH 是针铁矿（α-FeOOH）、四方纤铁矿（β-FeOOH）、纤铁矿（γ-FeOOH）的高压多形相。实验和数值模拟研究证实，ε-FeOOH可以在约10 GPa 条件下形成，低压下无法稳定存在[13-22]，与本计算结果（13 GPa）相符合。

Py-FeO2和Py-FeOOH 的弹性常数C11、C12、C44随压强的增加而线性增加（见图6），整体上看，Py-FeOOH 的弹性常数大于Py-FeO2的弹性常数。Py-FeO2和Py-FeOOH 的弹性常数变化与前人结果[9-10]基本一致。ε-FeOOH 的弹性常数随压强的变化如图6 和图7 所示，其中C11、C12、C22、C33、C66、C13、C23随着压强的增加线性增加，C44随压强的增加而减小，C55随压强的增加变化不大。ε-FeOOH 的弹性常数在33 GPa 附近大幅变化，这是由于结构随着压强变化而产生了氢键对称化现象。这一现象表明，ε-FeOOH 的弹性常数与氢键强度密切相关。对角弹性常数C44体现了晶格中（100）晶面沿着（010）晶向剪切应变的抵抗能力，（100）晶面上只有共边的氧原子和氢键，弹性常数C44在高压下的强度较低。压强的增加使C44逐渐减小，意味着额外的压力将导致ε-FeOOH 的结构稳定性降低。对于正交晶系结构的ε-FeOOH，需要满足的稳定性条件为：Cii＞ 0 （i= 1, 2, 3, 4, 5, 6），C11+C22+C33+ 2（C12+C13+C23） ＞ 0，C11+C22-2C12＞ 0，C11+C33- 2C13＞ 0，C22+C33- 2C33＞ 0。因此，ε-FeOOH 能稳定存在的压力边界（170 GPa）被限定。

图6 Py-FeO2、Py-FeOOH 和ε-FeOOH 的弹性常数（C11、C12、C44）随压强的变化Fig. 6 Elastic constants (C11, C12, C44) of Py-FeO2,Py-FeOOH and ε-FeOOH under high pressures

图7 ε-FeOOH 的弹性常数（C13、C23、C22、C33、C55、C66）随压强的变化Fig. 7 Elastic constants (C13, C23, C22, C33, C55, C66) of ε-FeOOH under high pressures

物质的弹性性质通常还可以用体积模量（K）和剪切模量（G）表示，一般采用Voigt 模型和Reuss 模型计算得到。由Voigt-Reuss-Hill 理论[50]可知，Voigt 模型和Reuss 模型计算得到的弹性模量为多晶体弹性模量的上、下限，而Hill 模型则是两者的数学平均，得到的计算结果与实验值的一致性更好，因此本研究采用Hill 模型计算弹性模量。

Py-FeO2、Py-FeOOH 和ε-FeOOH 的体积模量和剪切模量随压强的变化见图8。对比体积模量可知，Py-FeO2最高，Py-FeOOH 和ε-FeOOH 在高压下几乎一致；Py-FeO2的剪切模量最大，ε-FeOOH的剪切模量最小。Py-FeO2、Py-FeOOH 和ε-FeOOH的体积模量随着压强的增加线性增加，Py-FeO2和Py-FeOOH 的剪切模量随压强的增加线性增加。ε-FeOOH 的弹性常数C44随压强的增加而减小，剪切模量由弹性常数推导[45]可得（与C44的变化正相关），因而ε-FeOOH 的剪切模量呈现随压强增加而下降的趋势。ε-FeOOH 的弹性模量在33 GPa 附近大幅增加，这是由结构内的氢键对称化引起的。压力诱导ε-FeOOH 结构发生氢键对称化，提高了氢键的强度，致使弹性模量随之增加。Py-FeO2、Py-FeOOH 和ε-FeOOH 的体积模量随压强变化速率大于剪切模量随压强变化速率。

图8 Py-FeO2、Py-FeOOH 和ε-FeOOH 的弹性模量随压强的变化Fig. 8 Bulk and shear moduli of Py-FeO2, Py-FeOOH and ε-FeOOH under high pressures

2.3 波速随压强变化的关系

Py-FeO2、Py-FeOOH 和ε-FeOOH 是地幔及核幔边界的重要组成矿物，了解不同深度的地幔物质的物理性质，如密度、弹性性质和地震波速，可以更好地限定地幔物质的岩石学和化学组成模型。根据不同压强下Py-FeO2、Py-FeOOH 和ε-FeOOH 的密度和弹性模量计算得到的剪切波速（vs）和压缩波速（vp）如图9 所示。与实验数据[11,16]相比，Py-FeO2、Py-FeOOH 的剪切波速的相对偏差分别为32%和22%，压缩波速的相对偏差分别为15%和6.8%。NRIXS 实验与理论计算的绝对值相比有一定偏移，其原因在于在高压和高温条件下进行NRIXS 实验极具挑战性[11]。ε-FeOOH 的剪切波速的相对偏差约为6%，压缩波速与Zhuang 等[13]的研究结果基本一致。

图9 Py-FeO2、Py-FeOOH 和ε-FeOOH 的地震波速随深度的变化Fig. 9 Variations of seismic wave velocities of Py-FeO2, Py-FeOOH and ε-FeOOH with depth

Py-FeO2、Py-FeOOH 和ε-FeOOH 的压缩波速均随压强的增加而增加，Py-FeO2的剪切波速随压强的增加而增加，增加速率相对较小。Py-FeOOH 的剪切波速在0～2 000 km 深度范围内随深度的增加而增大，在2 000～6 000 km 深度范围内变化较小（5.8 km/s ＜vs＜ 6.0 km/s）；ε-FeOOH 的剪切波速在33 GPa（深度约900 km）发生突变，在0～900 km 深度范围内随深度的增加而减小，在900～3 000 km 深度范围内随深度的增加先增大后减小。ε-FeOOH 在33 GPa 压力下发生的氢键对称化现象可能在一定程度上影响地幔过渡带的密度和地震波速突变现象（从图9 中可以看到波速突变）。ε-FeOOH 的剪切波速在50 GPa 之后再次出现波速随压力增大而减小的情况，但其压缩波速却呈现随压力增大而增大的趋势，这是因为求解压缩波速时方程中的体积模量远远大于剪切模量，削弱了结构变化带来的影响。ε-FeOOH 的剪切波速发生了两次波速软化，可能与压力诱导下ε-FeOOH 的结构相变有关。Thompson等[28]的最新研究表明，在0～75 GPa 压力范围内，ε-FeOOH 有两次连续的结构相变，(18±1) GPa 时压力诱导的氢键对称化导致P21nm→Pnnm相变，在(45±2) GPa 经历等结构的一级相变。ε-FeOOH 的两次波速软化很可能是两次连续的结构相变导致的，具体分析还需要更多的数据支撑。

Py-FeO2、Py-FeOOH 和ε-FeOOH 的波速之间的大小关系显示，ε-FeOOH 的波速最低，而Py-FeO2的波速最高。在地球深部，水的存在对地震波速有较大影响，计算结果表明，不含水的Py-FeO2波速明显较高。Py-FeO2、Py-FeOOH 和ε-FeOOH 的波速曲线在400～700 km 深度范围内与PREM 和 AK135 波速模型相交，在700～2 898 km（核幔边界）深度范围内的波速远远低于 PREM 和 AK135 波速模型。

综合计算结果，可以发现，黄铁矿型FeO2和FeOOH 具有高密度、低波速的特点，与地幔ULVZs 的性质一致。地核-地幔边界区域的ULVZs 对地球的化学成分和热结构具有重要影响，但它们的起源一直存在争论[14]。如果Py-FeO2和Py-FeOOH 富集下沉到核幔边界，可能会成为ULVZs 的来源。

3 结 论

运用第一性原理方法，计算了ε-FeOOH、Py-FeOOH 和Py-FeO2在不同压强下的晶体结构变化和弹性性质。Py-FeO2和Py-FeOOH 的晶格常数随压强的增加呈现逐渐减小的趋势，晶胞体积也随着压强的增加而降低。ε-FeOOH 的晶格常数整体变化趋势是随着压强的增加逐渐减小，0～170 GPa 区间c轴长度相比a、b轴长度大幅减小，但在33 GPa 附近有异常突变，a、b轴长度突然增加，c轴长度突然减小。c轴上只有共边的氧原子和氢键，相比有共边的铁原子以及铁氧键的a、b轴，更容易被压力诱导压缩。Py-FeO2、Py-FeOOH 和ε-FeOOH 的晶胞密度随压强的增加而增加，Py-FeO2、Py-FeOOH 和ε-FeOOH 的密度在同等压强条件下依次减小，Py-FeO2的密度最高。

Py-FeO2和Py-FeOOH 的弹性常数C11、C12、C44随压强增加线性增加；ε-FeOOH 的弹性常数C11、C22、C33、C66、C12、C13、C23随压强的增加线性增加，C44随压强的增加线性减小，C55几乎不随压强变化。ε-FeOOH 的弹性常数在33 GPa 附近大幅变化，这是由于结构随着压强变化而产生的氢键对称化现象。这一现象表明，ε-FeOOH 的弹性常数与氢键强度密切相关。Py-FeO2、Py-FeOOH 和ε-FeOOH 的体积模量以及Py-FeO2、Py-FeOOH 的剪切模量随压强的增加线性增加，ε-FeOOH 的剪切模量在33 GPa 突变升高。对比体积模量可知，Py-FeO2最高，Py-FeOOH 和ε-FeOOH 在高压下几乎一致；Py-FeO2的剪切模量最大，ε-FeOOH 的剪切模量最小。

Py-FeO2、Py-FeOOH 和ε-FeOOH 的压缩波速随压强的增加而增加，Py-FeO2的剪切波速随压强的增加而增加。Py-FeOOH 的剪切波速在0～2 000 km 深度范围内随深度增加而减小，在2 000 ～6 000 km深度范围内变化较小（5.8 km/s ＜vs＜ 6.0 km/s）；ε-FeOOH 的剪切波速在33 GPa（深度约900 km）发生突变。Py-FeO2、Py-FeOOH 和ε-FeOOH 的波速对比显示，ε-FeOOH 的波速最低，Py-FeO2的波速最高。综合研究发现，黄铁矿型FeO2和FeOOH 具有高密度、低波速的特点，与地幔ULVZs 的性质一致。Py-FeO2和Py-FeOOH 在形成之后，可能会继续富集下沉到核幔边界，成为ULVZs 的来源。ε-FeOOH 在超过33 GPa 压强条件下发生的氢键对称化给ε-FeOOH 的结构带来显著变化，同时氢键对称化会影响原子间相互作用，进而改变ε-FeOOH 的弹性性质和地震波速。

本研究丰富了Py-FeO2、Py-FeOOH 相在高压下的弹性模量数据，为研究地球深部含铁的含水矿物提供了新的数据支持，同时，对于高温高压下Py-FeO2、Py-FeOOH 和ε-FeOOH 及其固溶体在下地幔的赋存形式还需要开展深入研究，以探索其在地幔深部的赋存和对地幔物质组成和结构的影响。