*欢迎大家去观看原文，并引用

标题：Molecular Dynamics Simulation Study on the SiO2−H2O System:New Insights on Water Dissolution Mechanisms in Melts

期刊：ACS Earth Space Chem.

网址：https://doi.org/10.1021/acsearthspacechem.4c00305

一、文章摘要

采用反应力场（ReaxFF）分子动力学（MD）模拟方法，模拟了在1473 K 和100 -700 MPa条件下硅熔体与纯水之间的两相平衡状态。研究了熔体中的水溶解度、含氢物质以及溶解 - 解析反应。主要发现包括：（1）熔体中的主要含氢物质为 −Si−O−H 和分子水；这与以往的观点一致。此外，我们发现随着压力的增加，水的含量增加速度比 −Si−O−H 更快。我们还观察到总含量低于1% 的许多次要含氢物质，如 H、O−H、−Si−O(−H)−Si−、−Si−O−H−O−Si−、−Si−OH2、−Si−O−H−O−H、H3O、H−O−H−O−H、−Si−O−H−OH2 和 (H2O)2。 （2）水的溶解过程分为三个步骤：水分子通过 H2O + −Si⇔ −Si−OH2与熔体接触，随后 −H2O基团通过 −Si−OH2⇔ −Si−O−H + H 分解，然后各种含氢物质主要通过 −Si−O−H + −Si⇔ −Si−O(−H)−Si−达到平衡。熔体中的水析出过程恰好是溶解过程的逆过程。 （3）计算得出的水溶解度从100 MPa时的0.82 wt%增加到700 MPa时的4.37 wt%。这些数值平均比实验值低2.27重量百分比（47%）[霍尔茨等人，《美国矿物学杂志》2000 年，第85卷，第682页]，但与最近使用经典分子动力学模拟对含水流纹岩熔体进行的计算结果一致[杜菲尔斯等人，《化学地质学》2020年，第533 期，第119300 号]。目前的研究结果揭示了水在这些熔体中溶解的反应路径，并有助于对成分更为复杂的岩浆以及含有各种挥发物（包括水、硫、氯、氟、硼和磷）的岩浆进行模拟，从而全面理解火山和岩浆-热液过程。

二、计算图文

图 1.（a）二氧化硅熔体 -水体系的初始配置。 （b）经过α形状构建后的配置（水相已省略，以便能清晰显示界面）。H为白色，O为红色，Si为黄色。

图 2.在1473 K 和400 MPa条件下，硅氧熔体 -水体系中（a）总能量和（b）水浓度随时间的变化情况，表明该体系在28纳秒后就达到了总能量和水浓度的平衡状态。此处，整个熔体在空间上被分为外层过渡层和内层熔体，详情见第3.1节。此外，在本文的修订阶段，该实验运行时间从33 纳秒延长至45纳秒，以重新检验其平衡状态，结果表明平衡状态依然良好（图 S9）。此外，我们还在6.5纳秒的状态点下进行了均匀的水合硅氧熔体实验，并观察到单相系统在约4纳秒时达到平衡（图 S10），这还得到了硅、氧和氢原子扩散行为的支持（图 S11）。

图 3.（a）在100 -700 兆帕压力下，氧 -氢键的径向分布函数；（b）在100 -700 兆帕压力下，硅 -氧键的径向分布函数。

图 4.熔体球体半径与压力之间的关系。

图 5.水在含水硅酸熔体中的溶解度。黑色的点代表霍尔茨等人测定的水溶解度。28蓝色和红色的点分别代表整个熔体和内部熔体的水溶解度。

图 6.水在含水二氧化硅熔体各层中的浓度。

表 1.在1473 开尔文温度及各种压力条件下，全二氧化硅熔体的水溶性模拟值与内二氧化硅熔体的水溶性模拟值之间的对比分析

图7.含水二氧化硅体系中水溶性随压力平方根的变化情况

图 8 。含氢的物质在含水的二氧化硅熔体中。黄色：硅；红色：氧；白色：氢。（a）−硅−氧−氢。（b）水。（c）氢。（d）氧−氢。（e）−硅−氧（−氢）−硅−。（f）−硅−氧−氢−氧−硅−。（g）−硅−羟基氢。（h）−硅−氧−氢−氧−氢。（i）氢−氧。（j）氢−氧−氢−氧−氢。（k）−硅−氧−氢−羟基氢。（l）（水）2。突出显示的球体表示含氢物质的组成原子。

图 9 。（a）结构水（−Si−O−H）与自由水（H2O）的溶解度与压力之间的关系。 （b）自由水与结构水的比例通常会随着压力的增加而增大。

表 2.含氢物质在二氧化硅熔体相中的反应途径

图 10.反应频率与压力之间的关系。

三、计算分析

本文采用Reactive Force Field分子动力学（ReaxFF MD）模拟方法，研究SiO₂熔体与水的两相平衡系统（1473 K，100–700 MPa），重点分析水在熔体中的溶解机制。通过大规模并行模拟（总计659 ns），揭示了水溶解的分子路径、溶解度行为以及熔体结构变化。模拟基于LAMMPS软件包实现，使用ReaxFF力场（Cu/Al/Si/O/H参数集），结合alpha-shape算法进行相分离分析，并通过几何准则识别H-bearing物种和反应路径。以下从科研角度详细解析模拟分析内容及计算细节。

3.1 模拟计算分析内容

本文通过ReaxFF MD模拟，系统探究了以下核心内容，结合热力学和动力学分析，提供分子层面的机制解释：

水溶解度与压力关系

模拟量化了水在SiO₂熔体中的溶解度随压力的变化。结果显示，溶解度从0.82 wt%（100 MPa）增至4.37 wt%（700 MPa），呈平方根依赖性（图7），符合Henry定律。

关键发现：溶解度低于实验值（Holtz et al., 2000），平均偏差47%，归因于力场对水密度的低估（图5）。这突显了模型优化的必要性。

实验基准数据引用

采用Holtz et al. (2000) 的实验结果作为基准：SiO₂熔体在1200–1350°C、101–600 MPa下，水溶解度范围为2.48–8.05 wt%（表1）。

补充材料S1（表S1）汇总了1931–2023年间的实验数据，涵盖不同成分熔体（玄武质至流纹质），证明本文压力-温度条件（1473 K≈1200°C）符合实验可比范围。

模拟结果与实验的定量对比

本文模拟值（0.82–4.37 wt%）系统性低于实验值（图5），平均偏差2.27 wt%（47%）。

归因分析：

ReaxFF力场对纯水密度预测偏低~3.32%（图S5），导致接触熔体的水分子数量不足。

与Dufils et al. (2020) 的经典MD结果一致（偏差~50%），共同指向力场对水-硅相互作用的描述不足。

H-bearing物种的分布与演化

识别了熔体中H-bearing物种的类型及比例：主要物种为结构水（-Si-O-H）和自由水分子（H₂O），占比 >99%；次要物种（<1%）包括H、O-H、-Si-O(-H)-Si-等（图8）。

压力影响：随压力增加，自由水比例上升更快（图9），表明高压促进分子水溶解，而结构水在 >300 MPa时增长放缓，暗示Si-O键合位点饱和。

物种组成的谱学支持

模拟识别主要物种为结构水（-Si-O-H）和分子水（H₂O），占比>99%（图8），与实验观测一致：

Stolper (1982) 红外光谱证实硅酸盐熔体中存在OH⁻和H₂O。

Shen & Keppler (1995) 原位高温高压光谱直接观测两类物种。

次要物种的实验间接证据

模拟发现<1%的次要物种（如H₃O、-Si-O(-H)-Si-），实验虽未直接检测，但Burnham (1975) 的水解反应模型（H₂O + Si-O-Si → 2Si-O-H）暗示中间态存在。

水溶解/解吸的分子机制

溶解过程分为三步（表2）：

接触：H₂O与Si接触形成-Si-OH₂（反应R2，频率40.7%）。

解离：-Si-OH₂解离为-Si-O-H和H（反应R7，频率0.11%），为速率控制步。

平衡：-Si-O-H与Si反应形成-Si-O(-H)-Si-（反应R1，频率57.76%）。

解吸为逆过程，高压下R1频率增长慢于R2（图10），表明自由水解吸更敏感。

熔体结构分层与过渡层效应

熔体相存在2.55 Å厚度的过渡层（L1-L3），水浓度显著高于内部熔体（L4-L5）（图6）。浓度梯度表明界面动力学不平衡，影响溶解度计算（仅内部熔体代表真实溶解度）。

熔体半径随压力增大（16.78 Å at 100 MPa → 17.17 Å at 700 MPa）（图4），归因于水溶解导致的体积膨胀。

3.2溶解机制的实验-模拟互证

三步溶解路径的合理性

传统实验模型（Burnham 1975, Stolper 1982）提出单步水解（反应式1），但模拟显示：

直接反应频率仅0.004%（图S13），而分步反应（R2+R7+R8）频率>40%（表2）。

实验无法捕捉瞬态中间体（如-Si-OH₂），模拟则揭示其作为关键过渡态（R2）。

压力效应的实验佐证

模拟显示分子水比例随压力上升（图9b），与Holtz实验数据中高压下分子水峰增强现象一致（IR光谱特征）。

3.3模型局限性与实验启示

力场优化的实验指导

模拟与实验的溶解度差异（图5）提示：需引入极化项（如AMOEBA）或优化ReaxFF的O-H键参数，以提升水密度预测精度（图S5）。

实验设计的模拟反馈

模拟发现熔体过渡层（图6）导致实验可能高估溶解度（若未区分界面/内部），建议后续实验采用微区光谱（如μ-FTIR）验证浓度梯度。

总结：计算与实验的协同价值

本文通过“模拟预测→实验对照→机制修正→模型优化” 的闭环，彰显分子模拟在补充实验盲区（如界面动力学、瞬态物种）中的不可替代性，为复杂挥发分（S, Cl）在岩浆中的行为研究提供新范式。

四、计算方法

4.1力场模型与参数

力场基础

采用ReaxFF力场（Cu/Al/Si/O/H参数集），源自Psofogiannakis et al. (2015)，继承自van Duin et al. (2003)的Si/O模型和Fogarty et al. (2010)的Si/O/H模型。

作用机制：通过动态键序计算化学键断裂/形成，可描述水分子解离、Si-O网络重组等反应（优于传统非反应力场）。

力场验证

密度与结构：对SiO₂熔体（3000 K, 1 bar）的预测密度（2.061 g/cm³）接近实验值（2.094 g/cm³），Si-O配位数（4.01）匹配FPMD数据（4.03）（图S1-S4）。

局限性：对纯水密度预测偏低~3.32%（图S5），扩散系数高估~38%（图S6），但对熔体-水界面反应适用（S2补充材料）。

4.2系统构建与初始化

初始配置生成

熔体制备：

从cristobalite晶胞（4 SiO₂）扩展至6×6×6超胞（864 SiO₂）。

在300 K、0.1 MPa下弛豫250 ps，后以8 K/ps速率升温至4300 K熔融1 ns。

切割500 Si/1000 O的球形熔体（半径~17.9 Å）。

两相系统：

熔体球置于45 Å立方盒中心，Packmol填充1000 H₂O分子（最小间距1.6 Å），避免原子重叠（图1a）。

边界与系综设置

周期性边界：三维周期性盒子。

系综：NPT（恒粒子数、压力、温度）。

控温/压方法：Nosé-Hoover恒温器（阻尼常数50 fs）和恒压器（阻尼常数25 fs）。

4.3 模拟流程与参数

平衡过程

温度：恒定1473 K（避免低温结晶或高温混相）。

压力阶梯：从100 MPa开始，以10 MPa/30 ps升至目标压力（100–700 MPa），后恒压运行至平衡（图2）。

平衡判定：系统能量与水浓度波动稳定 >2 ns（图2, S9）。

总时长：28组并行模拟，总计659 ns（平均23.5 ns/组）。

关键动力学参数

时间步长：0.25 fs（确保O-H键振动稳定性）。

轨迹输出：每50 fs保存一帧（用于物种与反应分析）。

4.4、数据分析方法

相分离与溶解度计算

Alpha-shape算法：

约束半径10 Å，识别熔体相为多面体（顶点为表面Si原子）（图1b）。

原子归属判定：若原子与质心连线不穿过多面体面，则为熔体内部原子。

溶解度公式：

S=0.5N_H*MH₂O/[N_H*M_H+N_O*M_O+(N_Si−N_Si^vertex)M_Si]×100%

其中N_Si^vertex为表面Si原子数（排除界面效应）。

H-bearing物种识别

键长阈值：

O-H键：≤1.35 Å（基于g_OH(r)第一谷位置，图3a）。

Si-O键：≤2.5 Å（覆盖g_SiO(r)第一谷内小峰，图3b）。

遍历算法：

以O原子为中心搜索键合原子，树状扩展至端点，分类物种（如-Si-O-H、H₂O、-Si-OH₂等）。

反应路径统计

帧间差异分析：

对比相邻轨迹帧（间隔50 fs），记录物种变化（如R1: -Si-O-H + Si ↔ -Si-O(-H)-Si-）。

频率计算：

累计最后2 ns的105,252次反应，按类型分类（表2）。

4.5关键参数总结表