白云母–流体界面盐离子分布的计算模拟研究: 对离子地球化学循环的启示

梁 箫, 陈 锰

白云母–流体界面盐离子分布的计算模拟研究: 对离子地球化学循环的启示

梁 箫1, 2, 陈 锰1*

(1. 中国科学院 广州地球化学研究所矿物学与成矿学重点实验室/广东省矿物物理与材料研究开发重点实验室, 广东 广州 510640; 2. 中国科学院大学, 北京 100049)

黏土矿物是表生土壤、沉积物、碎屑岩结晶颗粒表面的常见矿物, 常与地质流体处于共存状态。离子–水–黏土矿物的相互作用影响着表生环境的流体组成、元素输运、土壤有机质分馏和沉积物的力学特性等。通过分子动力学模拟认识离子在黏土矿物–流体界面的分布, 是揭示黏土矿物参与离子地球化学循环过程的基础。以白云母为研究对象, 使用分子动力学模拟和热力学计算方法, 分析了白云母–碱金属盐溶液界面的离子吸附结构和静电作用特性。研究发现Na+离子在白云母表面形成了两个内球络合吸附层, K+、Cs+离子仅形成一个内球络合吸附层。内球络合吸附层外可以形成多个外球络合吸附层, 表现为离子的层状吸附结构。增大盐浓度没有改变同种离子的吸附结构, 仅改变吸附层分布的大小。定量对比了使用伞状采样方法计算的离子吸附自由能曲线和基于离子密度分布得到的自由能曲线, 二者总体一致, 基本确证了模拟采样的充足性。研究表明界面流体不能视为均一的介电环境。

白云母; 流体; 碱金属离子; 自由能; 模拟

0 引 言

离子在矿物–流体界面的分配与自然过程息息相关。黏土矿物是表生土壤、沉积物、碎屑岩结晶颗粒表面的常见矿物, 常常与地质流体处于共存状态(Stoffelsma et al., 2010; Tinnacher et al., 2016; Resasco et al., 2017; Thomas and Peter, 2018; Magnussen and Gross, 2019; Waegele et al., 2019)。黏土矿物由带负电荷的片层组成, 并通过层间域或外表面的离子吸附实现电中和(Rotenberg et al., 2007; Lee et al., 2013)。黏土矿物表面的离子吸附和交换控制着含水层中的流体成分, 并影响有毒和放射性废物的地质储存(Rotenberg et al., 2007; Sakuma and Kawamura, 2011; Tinnacher et al., 2016; Waegele et al., 2019)。离子–水–黏土矿物的相互作用影响表生环境的流体组成、元素输运、土壤有机质分馏和沉积物的力学特性等(Bourg et al., 2017; Brown et al., 2017)。因此, 揭示黏土矿物的层间和外表面离子的竞争吸附、交换机理十分重要。过去通过结合分子动力学模拟和热力学计算, 揭示了黏土矿物层间域的离子交换机理(Zhou et al., 2020)。针对白云母的实验研究较多, 许多实验方法已经用于黏土矿物–电解质溶液界面离子的研究中, 如X射线反射、原子力显微镜等(Lee et al., 2013, 2017, 2021; Zou et al., 2021)。本研究以白云母为对象, 使用分子动力学模拟和热力学计算方法, 分析白云母–碱金属盐溶液界面的离子吸附结构和静电作用特性, 探讨黏土矿物外表面对碱金属离子的吸附和分配作用。

双电层结构(electric double layer, EDL)是常用来描述矿物表面离子行为的模型(Shchukin et al., 1958), 带正电的离子被吸附在黏土矿物表面形成吸附层(Stern层), 吸附层外为扩散层(Gouy-Chapman层)。根据Gouy-Chapman-Stern(GCS)模型, 静电势在吸附层内线性下降, 而在扩散层中呈指数下降。GCS模型依赖于点电荷假设和连续介电环境假设, 连续介电环境假设扩散层中的介电常数与纯溶液相同, 从实验数据和模拟结果可以看出GCS模型的局限性, 即使是相同电荷的离子, 不同的离子半径或离子类型会造成吸附层的结构差异(Strmcnik et al., 2009; Ricci et al., 2014; Lee et al., 2017; Resasco et al., 2017)。不同的碱金属离子, 虽然具有同样的离子电荷, 但在白云母–流体界面吸附时与水分子形成的络合结构存在差异(Bourg et al., 2017; Kobayashi et al., 2017)。经典密度泛函理论(density functional theory, DFT)可用于理论计算离子的密度分布和矿物–流体界面的静电势分布(Voukadinova and Gillespie, 2019; Cats et al., 2021; Wang et al., 2021), 这种方法优于传统的GCS模型, 因为它已经将离子尺寸的影响考虑在内。但是该理论仍然以均一的介电环境作为假设, 这与实际情况不符。

分子模拟被广泛应用于表界面的离子吸附性质的研究, 包括第一性原理分子动力学模拟和基于经验力场的分子动力学模拟。分子动力学模拟不以均一的介电环境为假设, 较DFT方法更优。使用第一性原理分子动力学方法模拟的固–液界面体系的尺寸比较小(Creazzo et al., 2018), 不能完全反映固液界面的性质, 因此, 基于经验力场的模拟方法被更广泛地应用于这一领域(Rotenberg et al., 2007; Sakuma and Kawamura, 2011; Dewan et al., 2014; Bourg et al., 2017; Kobayashi et al., 2017)。过去使用经验力场研究与黏土矿物表面结合的阳离子(包括Li+、Na+、K+、Cs+、Ca2+、Mg2+、Sr2+等), 这些离子在表面的吸附形态可以分为外球络合和内球络合。外球络合是指离子被水分子组成的水化壳包围, 然后间接吸附于矿物表面; 内球络合是指离子直接吸附于矿物表面, 离子与矿物表面间没有被水化壳隔断。在碱金属离子里, Li+、Na+离子倾向于形成外球络合, K+、Cs+离子则倾向于形成内球络合, 碱土金属族的Ca2+、Mg2+、Sr2+离子均倾向于形成外球络合。但是, 当黏土矿物的四面体结构层存在较多离子替代配位的时候, 由于表面电荷增加, 内球络合的离子会增多(Delville, 1991; Skipper et al., 1991; Chatterjee et al., 2000; Young and Smith, 2000; Hensen et al., 2001; Sutton and Sposito, 2001, 2002; Tambach et al., 2004; Whitley and Smith, 2004; Malikova et al., 2004, 2006; Liu and Lu, 2006; Mignon et al., 2010; Suter et al., 2012; Boek, 2014; Zhang et al., 2014; Dazas et al., 2014, 2015; Vao-soongnern et al., 2015; Yang et al., 2017; Rahromostaqim and Sahimi, 2018)。

本研究通过描述黏土矿物的ClayFF力场(Cygan et al., 2004)、离子的Horinek力场(Horinek et al., 2009)和水分子的SPC/E模型(Berendsen et al., 1987), 研究不同碱金属离子(Na+、K+、Cs+)在白云母/流体界面的分布, 使用伞状取样方法(Torrie and Valleau, 1977)计算离子从体相流体吸附到矿物表面的自由能, 将该自由能曲线与基于离子分布推导的自由能曲线进行定量对比, 并考虑不同离子浓度对结果的影响。通过模拟计算离子分布是广泛使用的方法, 但是, 由于模拟可能存在系综不遍历的问题, 模拟统计的离子分布可能是偏离实际的。将基于离子分布推导的吸附自由能和直接计算的吸附自由能做对比, 是确定系综是否遍历的重要手段, 也是本研究对今后模拟研究的方法学贡献。

1 模拟与计算方法

白云母的化学式为K+[Al2(Si3Al)O10(OH)2]–, 本次研究的白云母为2M1多型, 白云母中各个原子的坐标由Richardson使用X射线衍射确定(=0.52 nm,=0.90 nm,2.01 nm,90.0°,95.8°,=90.0°) (Richardson and Richardson, 1982)。根据Loewenstein规则, 四面体片的铝原子对硅原子的替换位不相邻(Loewenstein, 1954)。本研究构建了一个12×6的白云母超胞, 其在平面内尺寸为6.3 nm× 5.4 nm, 被置于一个长方体的模拟盒子中。3层云母片层堆叠于盒子的方向上, 我们将K+离子插入云母的层间域用于中和电荷。每一层云母的厚度大约为2.7 nm。模拟体系包含3层云母片层, 相当于方向上的1.5个超胞。由于2M1型白云母具有两个不同的(001)表面, 两个不同表面可能存在差异, 所以本研究构建了非整数个超胞来制造两个不同的表面。在方向上, 厚度约为10 nm的盐溶液层被放置于白云母片层上。模拟盒子的所有方向均使用周期性边界条件。白云母与盐溶液界面系统的示意图见图1。本研究模拟NaCl、KCl和CsCl 3种盐溶液, 用于比较不同浓度、不同类型的盐溶液的离子在白云母表面的吸附情况, 并使用ClayFF力场(Cygan et al., 2004)描述白云母片层中的原子, Horinek力场(Horinek et al., 2009)描述溶液中的离子, SPC/E模型(Berendsen et al., 1987)描述水分子。这些力场模型均使用勒让德·琼斯势, 在力场形式上是相容的。在过去的研究里, 使用这些力场组合模拟蒙脱石的结果与实验认识一致(Zhou et al., 2020)。

分子动力学模拟采用GROMACS 5.1软件包(Hess et al., 2008)。模拟的时间步长为1.0 fs, 模拟采用z系综: 温度保持在298 K, 使用Nosé-Hoover控温器(Nosé, 1983; Hoover, 1985)控制温度; 使用Parrinello-Rahman控压方法(Parrinello and Rahman, 1982), 控压过程中调整盒子在方向的尺寸,平面尺寸不变, 压力维持为一标准大气压。短程非键作用(范德华力和静电力)的截断半径为1.2 nm。使用particle-mesh Ewald(PME)方法(Essmann et al., 1995)描述长程静电作用。

将NaCl、KCl、CsCl 3种盐溶液设置0.5 mol/L、1.0 mol/L和1.5 mol/L 3个浓度梯度。通过置换溶液中的水分子, 将盐离子插入到溶液片层中, 并将一定数量的K+离子以相同的方式插入溶液来平衡白云母片层的负电荷。对每个系统进行10 ns的NPT系综模拟, 最后2 ns的模拟轨迹用来收集数据。对每个体系的平衡模拟重复5次。由于白云母/溶液体系具有2个界面, 结构分析基于10组数据的平均值。

Si、Al、Na、K、Cl、O、H原子分别由黄色、粉色、蓝色、紫色、青色、红色和白色的小球或点表示。

本次研究计算了阴阳离子从云母表面到溶液内部的平均力势(potential of mean force, PMF)。随机选择一个离子并将其沿法线方向从表面拉出(使用弹簧势, 劲度系数=500 kJ/mol), 速度为0.01 nm/ps。每20 ps记录一个结构, 共记录8个结构并用它们进行伞状采样模拟(Torrie and Valleau, 1977)。对每个伞状采样的窗口进行0.2 ns的随机分子动力学模拟, 使用弹簧势(=500 kJ/mol)使离子的位置围绕目标位置波动。上述操作在每个系统中重复了10次, 然后使用加权直方分析方法(weighted histogram analysis method, WHAM)分析数据并得到PMF曲线及其误差(Hub et al., 2010)。

2 结果与分析

2.1 离子密度分布

通过计算各种离子在模拟体系方向上的数量密度(图2、3), 可以清晰地反映不同盐浓度条件下, 离子在白云母表面流体中的分布情况。密度分布图的横坐标为离子的坐标位置与白云母表面位置之差, 其中白云母表面位置定义为白云母表面O原子的平均位置。阳离子的密度分布图(图2a～c)展示了与其他研究发现的位置相似的密度峰(Bourg et al., 2017; Kobayashi et al., 2017)。Bourg et al. (2017)和 Kobayashi et al. (2017)研究指出, 距离白云母表面0.3 nm之内的密度峰代表了阳离子和水分子、表面基团形成了内球络合物。阳离子密度分布反映了类似的内球络合结构。不同浓度体系中同种阳离子的分布形式在定性上是一致的(图3a～c), K+和Cs+离子在表面0.3 nm内只形成了一个显著的单密度峰; Na+离子则形成两个密度峰。Na+的第1和第2个峰分别代表了位于硅氧骨架表面复三方孔洞中的Na+和结构表面带电氧原子正上方的Na+(Bourg et al., 2017)。K+和Cs+离子的单密度峰均代表位于表面复三方孔洞中的阳离子。在浓度较高的情况下(如1.5 mol/L), Na+离子在距离表面约0.5 nm处形成了一个宽而小的峰(图3a), 这代表外球络合的离子水化结构。对相同浓度体系中不同阳离子的密度分布(图2a～c)进行比较时发现, 在0.5 mol/L、1.0 mol/L、1.5 mol/L浓度条件下, 峰值大小顺序均为: K+> Cs+> Na+; 浓度增大, Cs+与K+的峰值逐渐接近。随着盐浓度上升, Na+和K+的密度峰逐渐远离云母表面(图3a、b), 但在Cs+的密度分布图中观察不到这种现象(图3c)。

图2 离子在白云母–流体界面的密度分布

图3 离子在白云母–流体界面的密度分布

通过计算阴离子(Cl−)的密度分布情况可见(图2d～f、图3d～f), 在不同阳离子种类和盐浓度情况下, Cl−都会呈现出两个峰, 出现的位置大致相同。在不同浓度的同种盐体系中, 随着盐浓度的增加, 峰值逐渐增大(图3d～f)。相同浓度的不同阳离子盐体系中Cl−的分布情况见图2d～f, 两个Cl−峰出现的位置基本一样, 峰值的大小有如下关系: 第1个峰CsCl > KCl > NaCl, 第2个峰NaCl > KCl > CsCl。

2.2 离子脱附的平均力势

阳离子和阴离子从云母表面到流体内部的自由能变化, 记为+()和−(), 它们和密度分布的关系为:

式中:ρ()和−()分别为阳离子和阴离子在位置的密度;b +和b −为它们在体相的密度;为玻尔兹曼常数;为温度。利用2.1节得到的密度分布, 我们可以据公式(1)、(2)推算+()和−()。但是, 由于密度分布可能受模拟的随机性和采样不足的影响, 可能和实际情况存在偏差。为此, 同时使用了伞状采样的方法直接计算+()和−(), 并与以密度分布推算的结果进行对比。这样, 一方面验证密度分布结果的可靠性; 另一方面, 确定离子从表面到内部的自由能变化。

通过伞状取样得到的PMF曲线基本与基于密度分布推导得到的PMF曲线大致吻合(图4), 说明方法的可靠性。PMF曲线中的最小值点即表示离子移动到此位置时, 体系的能量最低、最稳定。图中的PMF最小值出现的位置与浓度基本无关, 而是取决于阳离子的种类, 这一发现与前文密度分布的结果一致, 可以说明不同的阳离子在云母表面有不同的吸附状态和偏好。PMF曲线的局部极小值点代表亚稳定的状态。从同种阳离子不同盐浓度的PMF曲线中可以看出: Na+的PMF曲线有4个极小值点, 而K+和Cs+则分别只有两个极小值点。通过对极小值点出现位置与密度分布的结果进行对比, 可以看出Na+的前两个极小值点分别与密度分布中的前两个峰相对应; K+与Cs+的第1个极小值点对应了密度分布中最明显的峰, 它们代表形成内球络合结构的阳离子。Na+的PMF曲线上的第3个极小值点、K+与Cs+的第2个极小值点对应于形成外球络合结构的阳离子。对于0.5 mol/L的NaCl的体系, Na+的PMF曲线还明显出现了第4个极小值点(图4a), 该极小值点随溶液的盐浓度增大而变得不明显(图4b、4c)。这说明了阳离子在白云母表面为多层吸附结构, 而PMF曲线较之密度分布图更明显地展示了这种多层吸附结构。Martin-Jimenez et al. (2016)研究中, 高分辨原子力显微镜也揭示过这种多层吸附结构。

图4 离子从白云母表面到溶液内部脱附的平均力势

此外, 通过计算还得到了对应体系中Cl−的PMF曲线, 发现Cl−在不同离子种类和浓度体系中都存在两个最小值点, 出现在=0.5 nm两侧, 这与密度分布的结果相符。Cl−离子的PMF曲线也展现了离子多层吸附的特点。

通过伞状取样得到的PMF曲线与基于密度分布推导得到的PMF曲线在定量上存在一定的偏差, 尤其是距离表面较近处的阳离子的内球络合吸附层(图4)。这种偏差可能是采样不足所导致的。如果增加伞状采样的重复次数和通过平衡模拟计算密度分布的次数, 这种偏差将可以被消除。

2.3 白云母表面流体的静电势

双电层理论和经典密度泛函理论对矿物–流体界面离子分布的描述皆以均一的介电常数()为假设(Shchukin et al., 1958; Voukadinova and Gillespie, 2019; Cats et al., 2021; Wang et al., 2021)。我们基于分子动力学模拟分析的结果来验证这个假设与实际情况的相符程度。

式中:0为真空中的介电常数;为元电荷量;z为原子或离子的带电荷量;ρ()为原子或离子的密度分布。

图5为静电势的分布图。图中可见,()随着表面距离增大, 先是到达峰值, 然后以近似指数方式下降。()则呈现波动式的衰减, 反映了水分子在白云母表面的分层结构, 是由离子水合形成的多层吸附所造成的结果。如果白云母表面的流体为均一的介电环境,()和()随距离的变化方式应当相似。二者在定性上的差别说明介电常数在界面上的分布为非均一的。本质上, 这种非均一的介电环境, 是由于离子与水分子结合, 在白云母外表面形成多层吸附。这种多层吸附的特性造成了介电环境的波动。

图5 静电势(f和ψ)在白云母–流体界面的分布

3 讨 论

在方法上, 本研究通过结合基于平衡模拟计算得到的密度分布和使用伞状采样方法得到的PMF曲线共同揭示了离子在白云母表面的分布状态。过去的模拟研究既有通过平衡模拟计算得到密度分布(Bourg et al., 2017), 也有通过采样方法计算自由能(Kobayashi et al., 2017), 但没有将二者的分析进行定量对比。本研究对比利用密度分布推导的PMF曲线与使用伞状采样方法得到的PMF曲线, 确证了在模拟采样充足的情况下, 二者皆可准确揭示离子在矿物表面的分布状态。PMF曲线较之密度分布图, 更好地展示了离子的外球络合吸附层的存在。这种分析方法的结合, 有望今后继续用于矿物/流体界面物质分布的分析。

本研究说明矿物–流体界面的介电常数并非是均一的, 这对双电层理论和经典密度泛函理论所应用的均一介电常数假设(Shchukin et al., 1958; Voukadinovaand Gillespie, 2019; Cats et al., 2021; Wang et al., 2021)提出了质疑。由于离子在矿物表面形成多层吸附结构, 每个吸附层的离子都与水分子在一定程度上形成了水化壳结构。这样, 将水环境视为均一的介电环境与实际情况不符。今后应当注意将介电性质随表面位置的变化纳入考虑中。当然, 仍然需要在理论上探索准确计算介电常数的方法。

通过实验方法可以认识矿物–流体界面离子的分布、离子与表面的结合状态等。分子动力学模拟的认识和实验认识形成对照, 这可以在现象与原理两方面确立矿物–流体界面的精细结构认识。本研究在分子动力学模拟基础上采用热力学计算, 可以用于计算不同离子在白云母表面的分配关系, 具有重要前景。

4 结 论

(1) 对于不同浓度盐体系Na+、K+、Cs+离子在白云母表面的吸附, Na+离子形成了两个内球络合层, 而K+、Cs+离子形成一个内球络合层, 内球络合层外可以形成多个外球络合层, 表现为离子的层状吸附结构。

(2) 增大盐浓度在定性上没有改变同种离子的吸附结构, 仅在定量上改变吸附层分布的大小。

(3) 白云母–流体界面的介电常数不能以均一的介电常数来表示, 介电常数应当是白云母表面距离的函数。

致谢:感谢西南科技大学张红平副教授和西北农林科技大学张弛副教授给出的专业评审意见, 对本文质量的提升有很大帮助。

Berendsen H, Grigera J R, Straatsma T P. 1987. The missing term in effective pair potentials., 91(24): 6269–6271.

Boek E S. 2014. Molecular dynamics simulations of interlayer structure and mobility in hydrated Li-, Na- and K-montmorillonite clays., 112(9–10): 1472–1483.

Bourg I C, Lee S S, Fenter P, Tournassat C. 2017. Stern layer structure and energetics at mica-water interfaces., 121(17): 9402–9412.

Brown K M, Poeppe D, Josh M, Sample J, Even E, Saffer D, Tobin H, Hirose T, Kulongoski J T, Toczko S. 2017. Theaction of water films at Å-scales in the Earth: Implications for the Nankai subduction system., 463: 266–276.

Cats P, Evans R, Hrtel A, Roij R V. 2021. Primitive model electrolytes in the near and far field: Decay lengths from DFT and simulations., 154(12), 124504.

Chatterjee A, Iwasaki T, Ebina T. 2000. A novel method to correlate layer charge and the catalytic activity of 2∶1 dioctahedral smectite clays in terms of binding the interlayer cation surrounded by monohydrate., 104(35): 8216–8223.

Creazzo F, Galimberti D R, Pezzotti S, Gaigeot M P. 2018. DFT-MD of the (110)-Co3O4cobalt oxide semiconductor in contact with liquid water, preliminary chemical and physical insights into the electrochemical environment., 150(4), 041721.

Cygan R T, Liang J J, Kalinichev A G. 2004. Molecular models of hydroxide, oxyhydroxide, and clay phases and the development of a general force field., 108(4): 1255–1266.

Dazas B, Ferrage E, Delville A, Lanson B. 2014. Interlayer structure model of tri-hydrated low-charge smectite by X-ray diffraction and Monte Carlo modeling in the Grand Canonical ensemble., 99(8–9): 1724–1735.

Dazas B, Lanson B, Delville A, Robert J L, Komarneni S, Michot L J, Ferrage E. 2015. Influence of tetrahedral layer charge on the organization of inter layer water and ions in synthetic Na-saturated smectites., 119(8): 4158–4172.

Delville A. 1991. Modeling the clay-water interface., 7(3): 547–555.

Dewan S, Carnevale V, Bankura A, Eftekhari-Bafrooei A, Fiorin G, Klein M L, Borguet E. 2014. Structure of water at charged interfaces: A molecular dynamics study., 30(27): 8056–8065.

Essmann U, Perera L, Berkowitz M L, Darden T, Lee H, Pedersen L G. 1995. A smooth particle mesh Ewald method., 103(19): 8577– 8593.

Hensen E J M, Tambach T J, Bliek A, Smit B. 2001. Adsorption isotherms of water in Li-, Na-, and K-montmorillonite by molecular simulation., 115(7): 3322–3329.

Hess B, Kutzner C, van der Spoel D, Lindahl E. 2008. GROMACS 4: Algorithms for highly efficient, load-balanced, and scalable molecular simulation., 4(3): 435–447.

Hoover W G. 1985. Canonical dynamics: Equilibrium phase- space distributions.,31(3): 1695–1697.

Horinek D, Herz A, Vrbka L, Sedlmeier F, Mamatkulov S I, Netz R R. 2009. Specific ion adsorption at the air/water interface: The role of hydrophobic solvation., 479(4–6): 173–183.

Hub J S, de Groot B L, van der Spoel D. 2010. g_wham—A free weighted histogram analysis implementation including robust error and autocorrelation estimates., 6(12): 3713–3720.

Kobayashi K, Liang Y F, Murata S, Matsuoka T, Takahashi S, Nishi N, Sakka T. 2017. Ion distribution and hydration structure in the stern layer on muscovite Surface., 33(15): 3892– 3899.

Lee S S, Fenter P, Nagy K L, Sturchio N C. 2013. Changes in adsorption free energy and speciation during competitiveadsorption between monovalent cations at the muscovite (001)-water interface., 123: 416–426.

Lee S S, Fenter P, Nagy K L, Sturchio N C. 2017. Real-time observation of cation exchange kinetics and dynamics at the muscovite-water interface., 8(1), 15826.

Lee S S, Koishi A, Bourg I C, Fenter P. 2021. Ion correlations drive charge overscreening and heterogeneous nucleation at solid-aqueous electrolyte interfaces., 118(32), e2105154118.

Liu X D, Lu X C. 2006. A thermodynamic understanding of clay-swelling inhibition by potassium ions., 118(38): 6448–6451.

Loewenstein W. 1954. The distribution of aluminum in the tetrahedra of silicates and aluminates., 39(1–2): 92–96.

Magnussen O M, Gross A. 2019. Toward an atomic-scale understanding of electrochemical interface structure and dynamics., 141(12): 4777–4790.

Malikova N, Cadene A, Marry V, Dubois E, Turq P. 2006. Diffusion of water in clays on the microscopic scale: Modeling and experiment., 110(7): 3206–3214.

Malikova N, Marry V, Dufreche J F, Simon C, Turq P, Giffaut E. 2004. Temperature effect in a montmorillonite clay at low hydration-microscopic simulation., 102(18): 1965–1977.

Martin-Jimenez D, Chacon E, Tarazona P, Garcia R. 2016. Atomically resolved three-dimensional structures of electrolyte aqueous solutions near a solid surface., 7(1), 12164.

Mignon P, Ugliengo P, Sodupe M, Hernandez E R. 2010.molecular dynamics study of the hydration of Li+, Na+and K+in a montmorillonite model. Influence of isomorphic substitution., 12(3): 688–697.

Nosé S C. 1983. A molecular dynamics method for simulations in the canonical ensemble., 52(2): 255–268.

Parrinello M, Rahman A. 1982. Strain fluctuations and elasticconstants., 76(5): 2662–2666.

Rahromostaqim M, Sahimi M. 2018. Molecular dynamics simulation of hydration and swelling of mixed-layer clays., 122(26): 14631– 14639.

Resasco J, Chen L D, Clark E, Tsai C, Hahn C, Jaramillo T F, Chan K, Bell A T. 2017. Promoter effects of alkali metal cations on the electrochemical reduction of carbon dioxide., 139(32): 11277–11287.

Ricci M, Spijker P, Voïtchovsky K. 2014. Water-induced correlation between single ions imaged at the solid-liquid interface., 5(1): 1–8.

Richardson S M, Richardson J W. 1982. Crystal structure of a pink muscovite from Archer’s Post, Kenya: Implications for reverse pleochroism in dioctahedral micas., 67(1–2): 69–75.

Rotenberg B, Marry V, Vuilleumier R, Malikova N, Simon C, Turq P. 2007. Water and ions in clays: Unraveling the interlayer/micropore exchange using molecular dynamics., 71(21): 5089–5101.

Sakuma H, Kawamura K. 2011. Structure and dynamics of water on Li+-, Na+-, K+-, Cs+-, H3O+-exchanged muscovite surfaces: A molecular dynamics study., 75(1): 63–81.

Shchukin E D, Pertsov A V, Amelina E A, Zelenev A S. Colloid and Surface Chemistry. 2001. Amsterdam: Elsevier: 165–220.

Skipper N T, Refson K, McConnell J D C. 1991. Computer simulation of interlayer water in 2∶1 clays., 94(11): 7434–7445.

Stoffelsma C, Rodriguez P, Garcia G, Garcia-Araez N, Strmcnik D, Markovic N M, Koper M T C. 2010. Promotion of the oxidation of carbon monoxide at stepped platinum single-crystal electrodes in alkaline media by lithium and beryllium cations., 132(45): 16127–16133.

Strmcnik D S, Kodama K, van der Vliet D, Greeley J P, Stamenkovic V R, Markovi N M. 2009. The role of non-covalent interactions in electrocatalytic fuel-cell reactions on platinum., 1(6): 466–472.

Suter J L, Sprik M, Boek E S. 2012. Free energies of absorptionof alkali ions onto beidellite and montmorillonite surfaces from constrained molecular dynamics simulations., 91: 109–119.

Sutton R, Sposito G. 2001. Molecular simulation of interlayer structure and dynamics in 12.4 angstrom Cs-smectite hydrates., 237(2): 174–184.

Sutton R, Sposito G. 2002. Animated molecular dynamics simulations of hydrated caesium-smectite interlayers., 3: 73–80.

Tambach T J, Hensen E J M, Smit B. 2004. Molecular simulations of swelling clay minerals., 108(23): 7586–7596.

Thomas G, Peter A E. 2018. Simulating donnan equilibria based on the Nernst-Planck equation., 232: 1–13.

Tinnacher R M, Holmboe M, Tournassat C, Bourg I C, Davis J A. 2016. Ion adsorption and diffusion in smectite: Molecular, pore, and continuum scale views., 177: 130–149.

Torrie G M, Valleau J P. 1977. Nonphysical sampling distributionsin Monte Carlo free-energy estimation: Umbrella sampling., 23(2): 187–199.

Vao-soongnern V, Pipatpanukul C, Horpibulsuk S. 2015. A combined X-ray absorption spectroscopy and molecular dynamic simulation to study the local structure potassium ion in hydrated montmorillonite., 50(21): 7126–7136.

Voukadinova A, Gillespie D. 2019. Energetics of counterion adsorption in the electrical double layer., 150(15), 154706.

Waegele M M, Gunathunge C M, Li J Y, Li X. 2019. How cations affect the electric double layer and the rates and selectivity of electrocatalytic processes., 151(16), 160902.

Wang X P, Liu K, Wu J Z. 2021. Demystifying the Stern layer at a metal-electrolyte interface: Local dielectric constant, specific ion adsorption, and partial charge transfer., 154(12), 124701.

Whitley H D, Smith D E. 2004. Free energy, energy, and entropy of swelling in Cs-, Na-, and Sr-montmorillonite clays., 120(11): 5387–5395.

Yang G M, Neretnieks I, Holmboe M. 2017. Atomistic simulations of cation hydration in sodium and calcium montmorillonite nanopores., 147(8), 084705.

Young D A, Smith D E. 2000. Simulations of clay mineral swelling and hydration: Dependence upon interlayer ion size and charge., 104(39): 9163–9170.

Zhang L H, Lu X C, Liu X D, Zhou J H, Zhou H Q. 2014. Hydration and mobility of interlayer ions of (Na, Ca)-montmorillonite: A molecular dynamics study., 118(51): 29811–29821.

Zhou H J, Chen M, Zhu R L, Lu X C, Zhu J X, He H P. 2020. Coupling between clay swelling/collapse and cationic partition., 285: 78–99.

Zou Y C, Mogg L, Clark N, Bacaksiz C, Milovanovic S, Sreepal V, Hao G P, Wang Y C, Hopkinson D G, Gorbachev R, Shaw S, Novoselov K S, Raveendran- Nair R, Peeters F M, Lozada-Hidalgo M, Haigh S J. 2021. Ion exchange in atomically thin clays and micas., 20(12): 1677–1682.

Molecular simulation study of distribution of salt ions at the muscovite-fluid interface: Implications to the ionic geochemical cycle

LIANG Xiao1, 2, CHEN Meng1*

(1. CAS Key Laboratory of Mineralogy and Metallogeny / Guangdong Provincial Key Laboratory of Mineral Physics and Materials, Guangzhou Institute of Geochemistry, Chinese Academy of Sciences, Guangzhou 510640, Guangdong, China; 2. University of Chinese Academy of Sciences, Beijing 100049, China)

Clay minerals are commonly found in soil and sediments. They generally coexist with geological fluid. The interaction between ions, water, and clay minerals affects fluid composition in the environment, transportation of elements, fractionation of organic matter, and mechanical properties of sediments. We investigated the distribution of ions at the clay mineral-fluid interface using molecular dynamics simulations to reveal the geochemical cycle of ions. We studied the ionic adsorption structure and electrostatic properties at the muscovite-alkali fluid interface through molecular dynamics simulations and thermodynamic calculations. We found that Na+ions formed two inner-sphere-complex adsorption layers on the surface of muscovite, whereas K+and Cs+ions only formed one inner-sphere-complex adsorption layer. The inner-sphere-complex adsorption layer was followed by several outer-sphere-complex adsorption layers, resulting in a layered adsorption structure of ions. Increasing the salt concentration did not change the adsorption structure when the salt type was maintained, even though the adsorption amount was altered. We compared desorption free energies derived using density distributions and the umbrella sampling method and found that they were consistent, which confirmed the sampling sufficiency. We found that the fluid on the surface cannot be regarded as a uniform dielectric environment.

muscovite; fluid; alkali ion; partition; simulation

P574; P579

A

0379-1726(2022)03-0261-10

10.19700/J.0379-1726.2022.03.002

2021-12-24;

2022-02-21

中国科学院青年创新促进会资助项目(2019345)资助。

梁箫(1996–), 男, 硕士研究生, 地质工程专业。E-mail: 740631451@qq.com

陈锰(1987–), 男, 副研究员, 主要从事计算矿物学研究。E-mail: chenmeng@gig.ac.cn