不同分子模型對甲烷水合物分解微觀特性表征

李佳，梁貞菊，王照亮，趙健，唐大偉

（1 中國石油大學（華東）新能源學院，山東青島266580； 2 大連理工大學能源與動力學院，遼寧大連116024）

引 言

甲烷水合物是在低溫高壓下形成的一種類似于冰的晶體包合物[1-2]，它以固態形式賦存于深海沉積物和永久凍土地區[3-4]。水分子通過氫鍵形成晶格主體，甲烷分子被包裹在水分子形成的籠形結構中，甲烷分子與水分子之間通過范德華力作用結合[5-6]。據估計，全球儲存在海底和永久凍土水合物中的甲烷約為1×1015～1.2×1017m3，全球天然氣水合物的潛在產量非?？捎^。由于其特殊結構，天然氣水合物對溫度壓力變化非常敏感，極易分解[7]。目前，國內外尚未實現甲烷水合物的大規模開發，主要原因是還未實現對甲烷水合物的可控分解。因此，研究水合物的分解動力學特性具有重要意義，可為水合物的開發利用提供重要依據。

近年來，國內外對水合物的分解過程進行了大量的實驗和理論研究[8-13]。研究發現，水合物的分解反應涉及氣、液、固三相變化，是一個復雜的吸熱過程。在實驗上很難從宏觀角度測定甲烷水合物分解的反應機理，而分子動力學模擬能夠克服實驗的局限性，可為水合物分解機理的微觀研究提供強有力的支持。國內外學者已采用多種主、客體分子模型進行了水合物分解的分子動力學模擬研究[14-15]。常見的應用于甲烷水合物分解的客體分子模型主要有OPLS-UA[16]聯合原子模型及TKM-AA[17-18]全原子模型。常見的水分子模型主要有TIP3P、SPC/E、TIP4P、TIP4P/Ice、TIP4P/2005、TIP4P-ew 和TIP5P 模型。由于TIP3P 模型很難正確反映水合物的結構特性，而TIP5P 模型存在計算量過大等缺點，目前在水合物的模擬中很少采用TIP3P 和TIP5P 模型[17]，廣泛應用的主體分子模型主要為SPC/E 模型及四點水分子模型。然而，由于力場參數不同，不同水模型預測的水的物理性質也會存在很大差異。因此，了解哪些水模型適用于在給定的熱力學條件下測量某一特性非常重要。對于某一特定的水模型，其力場參數不會隨著熱力學條件的變化而變化。因而單一的水模型無法同時適用于描述水合物生成(低溫)以及分解(高溫)條件下的兩種不同的熱力學狀態[19]。因此，需要兩種不同的水模型分別模擬水合物的生成和分解兩個獨立的過程。根據文獻可知，TIP4P/Ice[20-21]、TIP4P/2005[21-23]和TIP4P-Ew[24]水模型是天然氣水合物生成研究的首選模型。而對于水合物的分解特性研究多采用TIP4P[19,25]和SPC/E[16,18,26-28]模型研究水合物的分解動力學。即便如此，不同模型預測的水合物的屬性依然存在很大差異。因此，尋找適合于天然氣水合物分解研究的主客體分子模型，通過使用不同的水模型對水合物分解進行比較研究，并確定水合物分解的最佳水模型，仍然是一個需要關注的模擬研究領域。

綜上所述，本文分別采用水分子-甲烷的TIP4P-UA、TIP4P-AA、SPC/E-AA 和SPC/E-UA 模型模擬甲烷水合物的分解動力學特性。首先，模擬水合物分解過程中平衡溫度、序參數、不同區域勢能、質量密度以及逸出水合物的甲烷分子數目的變化等分解動力學參數，然后對水、界面、水合物區域進行了劃分，計算出分解過程的活化能，進而對比不同水分子、甲烷分子模型對分解的影響并研究不同模型下分解的微觀動力學機制。

1 計算方法

1.1 模擬體系

首先，構建SI 型甲烷水合物晶胞，其空間群為Pm3n，晶格長度為11.877 ?(1 ?=0.1 nm)，根據X 射線衍射實驗的結果確定C、O 原子的空間位置[29]，利用Bernal-Fowler 法則為C、O 原子添加H 原子以滿足整體凈偶極矩為零[30]。然后，將甲烷水合物晶胞沿X、Y、Z 方向擴展為3×3×6 的SI 型甲烷水合物超晶胞結構，沿Z 軸在甲烷水合物的兩側各放置1513個液態水分子。構建的XZ 平面內水-甲烷水合物-水初始模擬結構如圖1 所示，X、Y、Z 方向均采用周期性邊界條件, ABC 區域將在后文2.6 節詳細介紹。

圖1 初始模型Fig.1 Initial configuration of simulation system

1.2 模型選取

對水分子，選擇SPC/E、TIP4P 兩種水分子模型進行計算。其中，SPC/E 是三點模型，三個作用位置分別對應水分子的三個原子，三點模型具有計算簡便高效的特點。TIP4P 為四點模型，該模型假定有一個虛擬原子位于H-O-H 角平分線上，距離氧原子0.125 ?。四點模型更為復雜，可以更合理地描述液態水的結構和熱力學性質[31]。由于SPC/E、TIP4P模型均為剛性模型，模擬中需固定水分子的鍵長和鍵角。甲烷分子作為客體分子，最為簡單的模型為單一作用點的OPLS-UA聯合原子模型。此外，對甲烷分子的模擬還選用了TKM-AA 全原子模型，該模型有五個相互作用點，分別對應甲烷分子的五個原子。各種模型的參數如表1所示。

表1 水和甲烷模型參數Table 1 Parameters of water and methane models

分子間相互作用勢包括Lennard-Jones(L-J)勢及靜電勢，如式(1)所示

不同原子間的L-J 勢參數由Lorentz-Berthelot混合規則確定，如式(2)所示

1.3 模擬方法

采用LAMMPS 軟件進行模擬。采用速度-Verlet 算法積分牛頓運動方程，設置時間步長為1 fs，分子間相互作用的截斷半徑設置為12 ?，采用PPPM(particle-particle particle-mesh)法處理長程庫侖力，計算精度為10-6。系統的初始溫度設置為300 K，壓力設置為10 MPa，模擬步驟如下。

(1)對水合物做固定處理，將液態水置于NVT 系綜，采用Nose-Hoover 熱浴使之達到目標溫度，溫度的阻尼系數設置為0.2 ps，運行200 ps；

(2)繼續固定水合物，將液態水轉入NPT 系綜，設置目標壓力為10 MPa，繼續弛豫使之達到目標溫度，溫度的阻尼系數設置為1 ps，壓力的阻尼系數設置為2.5 ps，運行200 ps；

(3)解除對水合物的固定，將整個系統置于NPT系綜進行控溫控壓處理，運行30 ps；

(4)將整個系統置于NVE系綜分解，運行2 ns。

2 結果與討論

2.1 甲烷水合物的平衡溫度及分解熱

隨著分解的進行，水合物不斷從液態水吸收熱量，液態水的動能逐漸降低，溫度隨之下降。圖2為分解過程中溫度隨時間的變化。由圖2 可見，分解初始階段可用于提供熱量的動能比較充足，液態水溫度下降較快，隨著可用動能的消耗，液態水可提供的熱量有限，溫度緩慢下降。類似的溫度變化趨勢在Alavi 等[18]的模擬中也得到了驗證。此外，在恒定體積條件下，甲烷分子從籠形結構中逸出導致體系壓力升高，這也可能導致水合物進入相穩定區，從而使得分解過程終止。達到穩定狀態時，液態水的溫度稱為平衡溫度。此處采用式(3)計算體系的平衡溫度[25]。

針對圖2 中的溫度分布(黑色曲線)，每20 個數據點取其平均值，作出溫度隨時間的波動曲線(紅色曲線)，然后通過式（3）進行擬合，可以得到體系的平衡溫度(藍色虛線)。由圖2 可見，在初始溫度為300 K、初 始 壓 力 為10 MPa 的 條 件 下，TIP4P-AA 和TIP4P-UA 的平衡溫度較高，分別為(292.63±3.71)K和(291.86±1.55)K。SPC/E-AA 和SPC/E-UA 預測的平衡溫度分別為(286.36±4.33)K 和(282.89±3.62)K。由于設置的模擬初始溫度較高，所以在相同的壓力下，平衡溫度越高，水合物越容易達到平衡狀態不再分解，因而平衡溫度越高，結構越穩定[32]。Goel[33]測得10 MPa 甲烷水合物的平衡溫度為286.11 K，可見當水分子采用SPC/E 模型時得到的平衡溫度與實驗值更加接近。相比SPC/E-UA，SPC/E-AA 模型得到的平衡溫度更精確，這是因為TKM-AA 全原子模型對CH4分子的各個原子均設置點電荷，以再現CH4氣體分子的電八極矩，同時將L-J 勢分配給CH4的中心碳原子，而對氫原子不設置L-J相互作用，與OPLS-UA 模型相比，更能精確地對甲烷分子進行描述。

由于甲烷水合物的分解為吸熱反應，其反應方程式可以表示為

其中甲烷水合物的分解熱ΔHd可由Clausius-Clapeyron方程[34]求得，如下所示

式 中，Z 為 壓 縮 因 子，可 由Peng-Robinson 方程[35]求得，結合p-T 的對應關系[36]，進而可以求得采用TIP4P-AA、TIP4P-UA、SPC/E-AA 和SPC/E-UA模型時甲烷水合物的分解熱分別為53.85、53.27l、54.79 及56.68 kJ/mol。Rueff 等[37]運用DSC 量熱法測得的甲烷水合物的分解熱為54.49 kJ/mol，孫志高等[38]由恒溫壓力搜索法測得的水合物的分解熱為56.22 kJ/mol，由此可見當水分子采用SPC/E 模型計算甲烷水合物的分解熱時結果與實驗值更加吻合。

2.2 分解過程F3序參數的變化

圖2 不同模型的平衡溫度Fig.2 Equilibrium temperature under different models

在模擬過程中，甲烷水合物逐層向內分解，分解前沿區域由外向內推移，水合物-液態水界面向水合物側移動，可以通過計算氧原子的序參數[26]來區分類水合物水和液態水，進而確定水合物-液態水的界面位置。序參數F3的定義如下

由于水合物相水分子通過氫鍵形成四面體式排列，在理想水合物結構中鄰近水分子的氧原子形成的夾角為104.25°。因此，可以通過任一時刻鄰近水分子中氧原子夾角的測量，來衡量水合物中氫鍵網絡偏離理想四面體式排列的程度。理想狀態下水合物相水分子的F3約為0，而液相中水分子的排列具有無序性，F3約為0.1[39]。圖3 為不同模型下體系中水分子的F3沿Z 軸的分布情況?？梢耘卸‵3≈0.1 的區域為液態水區域，F3≈0 的區域為水合物區域，在液態水區域至水合物區域的過渡段F3出現跳躍式下降，該過渡段被界定為界面區域。圖3 中用虛線對液態水、水合物及界面區域做了區分?？梢钥闯?，在初始時刻t=0 時水合物區域F3接近于0，水合物兩側界面區域分別位于35.3～41.5 ? 和到109.4～116.9 ?。分解至t=400 ps 時，采用TIP4PAA 模型組合時界面區域分別位于40.3～47.1 ? 和98.7～105.2 ?。采用TIP4P-UA 模型組合時界面區域分別位于41.83～46.9 ? 和101.6～106.7 ?，同時水合物區域的F3較初始時刻有所上升，這是由于隨著分解的進行，液態水區域向水合物區域傳遞熱量，水合物內部水分子活躍程度增加，籠形結構出現扭曲變形造成的；而采用SPC/E-AA 模型組合時界面區域分別為48.7～55.1 ? 和88.7～94.8 ?，與TIP4P-AA 模型相比，該模型組合下，界面區域向內推移，且水合物區域F3與TIP4P-AA 模型相比略微升高。這主要是由于TIP4P 水模型與SPC/E 水模型所攜帶的電荷量不同引起的。TIP4P 水中氧原子與氫原子的電荷量大于SPC/E 水中氧原子與氫原子的電荷量，氧原子與鄰近氫原子之間的庫侖力更大，形成的氫鍵更穩固，因而籠形結構不易遭到破壞；此外，采用SPC/E-UA 模型組合時界面區域分別為57.1～62.1 ? 和97.3～103.6 ?。Sloan 等[1]預測的水合物-水的界面厚度約為5 ?，而Myshakin 等[40]估算的界面厚度在5～7 ? 之間，以上通過F3預測的各水合物-水界面區域厚度均在5～7.5 ? 之間，這與Sloan 等[1]及Myshakin 等[40]的結果相吻合，從而也驗證了此方法的正確性。

圖3 氧原子序參數沿Z軸變化Fig.3 F3 distribution along Z axis

圖4為不同模型下水合物中水的序參數F3在分解過程中的變化特征。由圖4 可見，四種不同模型下甲烷水合物在分解達到2000 ps 時整體的F3均未達到0.1，說明至模擬結束時依然有部分水合物尚未分解。在0～2000 ps 內，TIP4P 模型預測的F3變化比較平緩，分解結束時增至0.04 左右，說明采用該模型時水合物中水分子的有序度更高，甲烷水合物分解相對較少。采用SPC/E模型時氧原子的F3曲線斜率大于采用TIP4P 模型時氧原子的F3曲線斜率，說明SPC/E 模型下甲烷水合物的分解速率大于TIP4P 模型的預測結果?？蓮膭莺瘮到嵌确治鯯PC/E 模型的快速分解機理。由表1 可知，SPC/E 模型的σ與TIP4P模型相差很小，然而與TIP4P相比，SPC/E模型的勢阱深度ε 更小，因而分子間范德華力作用較小，同時氧原子與氫原子攜帶電荷較少導致氫鍵作用小，因而SPC/E 模型中單個水分子受到的束縛更小，水分子更加靈活。保持水分子模型不變，對比甲烷分子的AA 及UA 模型，發現F3變化特性基本一致，說明對甲烷分子的模擬選用TKM-AA 模型或OPLS-UA 模型對分解速率影響不大，這也表明在水合物的分解過程中水-水相互作用占主導地位。此外，圖4 還給出了水合物中水分子/甲烷分子的協同作用對分解的影響，由圖4可知，保持水分子模型及甲烷分子模型不變，當水分子/甲烷分子協同作用程度增大為原來的1.5 倍時，F3曲線斜率明顯降低，因而水合物分解更慢。這主要是因為當主客體分子協同作用增強時，客體分子的局域化特征增強，主客體分子之間的耦合程度增強，能量傳遞阻力增大，因而籠形結構穩定性增強。

圖4 水合物中F3序參數隨時間的變化Fig.4 Time dependence of F3 in methane hdyrate

2.3 甲烷分子勢能的變化

由于每個甲烷分子受到的勢能源于模擬盒中所有其他分子對它的相互作用，根據圖3 中液態水區、水合物區、界面區的界定標準，分別對分解過程三個不同區域的單個甲烷分子的平均勢能進行了計算。圖5 所示為單個甲烷分子勢能隨時間的變化。由圖5可見，隨著分解的進行，甲烷分子的勢能逐漸增加，這是由于NVE 分解過程體系的動能降低轉化為水合物分解所需熱量，從而導致了勢能的增加。對比不同區域的勢能可知液態水區單個甲烷分子勢能最大，界面區單個甲烷分子勢能次之，水合物區單個甲烷分子勢能最??；同時，不同區域甲烷分子勢能的增量也不一致，這也體現了溫度分布的不均勻性，即分解過程中不同區域存在溫度梯度和熱流梯度，從而驗證了水合物分解過程傳熱的非平衡特性[18]。對比不同模型預測的勢能可知，SPC/E-AA 模型預測的單個甲烷分子勢能大于TIP4PAA 模型單個甲烷分子勢能。甲烷分子勢能越低，穩定性越高，因而TIP4P分子模型穩定性越高，這也導致了圖4中TIP4P模型預測的分解速度較慢。

2.4 分解過程密度分布的變化

圖6為初始時刻及分解結束時不同分子模型的構象以及水合物中甲烷、水沿Z 軸的密度分布。在初始時刻，水合物中甲烷和水的密度分布比較規律，水的密度按照單雙峰呈現周期性變化，甲烷的密度按照高低峰呈現周期性變化。由圖6(b)ⅰ可知，對于完整的3×3×6型SI甲烷水合物而言，甲烷的密度分布共有12對高低峰組合。由于尚未分解，在左右兩側液態水區域甲烷和類水合物水的密度均為0。隨著分解的進行，籠形結構不斷破裂，甲烷分子、水分子開始從水合物中逸出，并擴散到液態水中，甲烷和水分子密度分布的規律性逐漸減弱。由圖6(b)ⅱ～ⅴ可見，隨著分解的進行，液態水區甲烷和水的密度開始出現非零正值。分解結束時，對于TIP4P-AA 模型和TIP4P-UA 模型，由圖6（a）Ⅱ～Ⅴ可知，尚有四層水合物籠形結構未分解，而對于SPC/E-AA 和SPC/E-UA 模型則只剩下約兩層沒有分解。因而，圖6(b)ⅱ～ⅴ中甲烷密度的高低峰殘存數分別為8 對、8 對、4 對、4 對。從圖6(a)Ⅳ所示的SPC/E-AA 模型預測的2 ns 時的構象圖可知，甲烷分子發生聚集形成甲烷氣泡，這與圖6(b)ⅳ中甲烷的密度在34? 和108? 附近呈現峰狀結構而水的密度在此處呈現谷狀結構相對應。這種現象的形成是由于隨著分解的進行大量的甲烷分子快速釋放到液態水中，當甲烷濃度超過液態水中甲烷的溶解度時，甲烷分子開始聚集形成甲烷富集區。

圖5 單個甲烷分子勢能隨時間的變化Fig.5 Variation of potential energy of single methane molecule with time

2.5 分解過程甲烷分子數的變化特征

由于甲烷水合物中十四面體籠形結構的直徑約為8 ?，結合Alavi 等[18]的研究方法，如果甲烷分子任一時刻的坐標與初始時刻的坐標距離超過8?，則認為甲烷分子已經從破碎的籠形結構中逃逸出來。圖7為分解過程逸出的甲烷分子數的變化特征。由圖可見，SPC/E-AA 模型、SPC/E-UA 模型預測的逸出甲烷分子數大于TIP4P-AA 模型及TIP4P-UA 模型的預測結果，SPC/E-AA 模型與SPC/E-UA 模型預測的逸出甲烷分子數目基本保持一致，這意味著采用SPC/E 模型時水合物的分解速度快于TIP4P 模型，而甲烷分子模型對分解的影響甚小，這與根據圖4得出的結論一致。從另外三種模型下逸出甲烷分子數的變化曲線可以看出，分解初始階段甲烷的逸出速率較大，之后逐漸減緩。這是因為該模擬在NVE 絕熱系綜進行，水合物分解是一個吸熱過程[25,41]，隨著分解的進行，NVE 系綜溫度逐漸下降，水合物分解所需驅動力不足，傳熱傳質阻力增大，因而分解速率不斷衰減。分解結束時，SPC/E-AA和SPC/E-UA 模型預測的逸出甲烷分子數約為290個，分解比例約為2/3，這與圖6 中分解結束時兩種模型的預測結果吻合較好，其中甲烷密度的殘存峰值對數為4對(殘存比例1/3)。對于TIP4P-AA 模型，紅色虛線標記的位置處逸出的甲烷分子數出現下滑，這主要是由于甲烷水合物在分解過程中發生了二次生成效應，當甲烷分子濃度高且存在過冷驅動力時逸出的甲烷分子會重新與水分子結合再次生成甲烷水合物。

根據Arrhenius 公式，甲烷水合物的分解速率可由甲烷分子數隨時間的變化來表示，如式(7)所示

由質量守恒方程可得甲烷分子數與分解界面之間的關系為

將式(8)代入式(7)可得

根據式(7)～式(9)，計算可得甲烷水合物分解所需的活化能，其中采用SPC/E-AA 模型預測的活化能為83.8 kJ/mol，采用SPC/E-UA模型預測的活化能為84.1 kJ/mol，采用TIP4P-AA 模型預測的活化能為86.7 kJ/mol，采用TIP4P-UA 模型時活化能為86.2 kJ/mol。Clarke 等[42]測得甲烷水合物分解所需的活化能為81 kJ/mol，Myshakin 等[40]通過對水合物在不同溫度下反應速度的擬合求得活化能為(82.2±2.1)kJ/mol，因而在NVE 系綜內采用SPC/E 模型時所預測的活化能與上述測量值吻合較好。

2.6 擴散系數

甲烷水合物在NVE 系綜分解時，原子從初始位置開始不停移動，ri(t)表示t 時刻時粒子的位置，粒子位移平方的平均值稱為均方位移，即

當t 值非常小時，R(t)呈指數增加；隨著t 值的增大R(t)近似直線分布。根據Einstein的擴散定律

由式(11)可知，通過均方位移的斜率可求出擴散系數。通過對甲烷水合物相平衡溫度及分解所需活化能的計算可知，當水分子采用SPC/E 模型時的計算結果與實驗值更加吻合，而TKM-AA 及OPLS-UA 模型對分解影響很小，因而此處采用SPC/E-AA 模型研究不同位置甲烷的擴散規律。針對SPC/E-AA 模型預測的水合物的分解過程，將完整的3×3×6 水合物分為如圖1 所示的A、B、C 三個區域，A 區域為最外層水合物，C 區域為最內層水合物，B 區域位于A 區域和C 區域之間。對這三個區域內甲烷分子的均方位移進行計算，可得到圖8 所示的不同區域甲烷的均方位移隨時間的變化關系。由圖8 可見，不同區域甲烷的均方位移存在很大差別，說明分解過程中不同區域存在濃度梯度，體現了分解過程傳質的非平衡特性[18]。與水接觸最近的A 區域內甲烷分子的均方位移最大，對應的甲烷分子擴散系數為1.4244×10-5cm2/s。與A 區域相比，B區域甲烷分子的擴散系數小一個量級，僅為6.6487×10-6cm2/s，因而B 層的傳質阻力要遠大于A 層。這主要是因為隨著A 層的分解，在水合物和水的界面處形成一層準液態膜，阻止B層內甲烷分子的擴散。此外，由圖6(a)Ⅳ體系構象圖中可見，分解結束時C區域尚未分解，因此甲烷分子只是在平衡位置附近振動及轉動，其均方位移基本恒定在一個稍微大于零的值附近，均方位移隨時間的變化曲線斜率為0，這表明C 區域甲烷的擴散系數為0。此外，由菲克定律可知：M = -D?ρ。其中等號左側項代表了質量的傳遞率，右側第一項代表擴散系數，右側第二項代表質量濃度梯度，即傳遞過程的推動力。從圖6 中甲烷的質量密度的變化規律可以看出A 區域甲烷的質量密度梯度最大，B 區域次之，因而甲烷水合物分解過程中B 區域分解時的推動力小于A 區域的分解推動力，同時根據圖8 的計算結果，A區域擴散系數最大，B 區域次之，進而導致A 區域分解時質量的傳遞率最大，B 區域的質量傳遞率小于A 區域。這也驗證了圖7 逸出的甲烷分子隨時間的變化曲線中起始段斜率較大，之后逐漸減小的變化趨勢。

圖6 體系構象和質量密度Fig.6 System configuration and density distribution

圖7 逸出籠子的甲烷分子總數隨時間的變化Fig.7 Number of methane molecules fleeing out of methane hydrate as function of time

圖8 SPC/E-AA模型不同區域甲烷均方位移隨時間的變化Fig.8 Mean square displacement of methane at different zone under SPC/E-AA model

3 結 論

(1) 在初始溫度300 K、初始壓力10 MPa 條件下，采用SPC/E 模型模擬得到的NVE 系綜內水合物分解的平衡溫度、活化能、分解熱與實驗測量結果更加接近。

(2)根據不同相態水的F3序參數，給出了區分液態水、水合物水及界面的方法，判定的液態水-水合物界面區域厚度為5～7.5 ?。對液態水、界面、水合物區甲烷的勢能進行計算，發現甲烷分子在水合物區的勢能最低，界面區勢能次之，液態水區勢能最高，體現了NVE系綜分解過程的非平衡傳熱特征。

(3)通過采用不同模型對分解過程中F3序參數、甲烷和水的密度分布及甲烷的逸出規律的模擬對比，說明甲烷分子模型選用TKM-AA 或OPLS-UA模型對分解的影響甚微，而不同水分子模型對分解影響較為明顯。由于TIP4P 模型下分子勢能較低，結構較穩定，相同溫度、壓力條件下SPC/E 模型預測的分解速率大于TIP4P模型。

(4)計算得到不同區域甲烷分子的擴散系數范圍為10-6～10-5cm2/s，不同區域甲烷的擴散系數差距明顯，由外向內傳質阻力逐漸增加，體現了NVE 系綜下分解過程的非平衡傳質特征。

符 號 說 明

As——甲烷水合物與液態水接觸的截面積，m2

cv——水的比定容熱容，J/(kg·K)

D——擴散系數，cm2/s

ΔHd——水合物的分解熱，kJ/mol

K0——反應常數

I——分解界面的位置，?

M——質量通量密度，kg/(m2·s)

MSD——均方位移，?2

mw0——液態水的初始質量，kg

nCH4——甲烷分子數

p——平衡壓力，Pa

Q+,Q-——原子所帶的正負電荷量，e

qi,qj——原子所帶電荷，e

R——摩爾氣體常數，J/(mol·K)

ri(t)——t時刻時粒子的位置

T——平衡壓力對應的平衡溫度，K

Teq——體系的平衡溫度，K

t——時間，ps

U(rij)——原子i與j之間的相互作用勢，4.18×103J/mol

Z——壓縮因子

ε——勢阱深度，4.18×103J/mol

ε0——真空介電常數，F/m

θjik——原子j，i，k形成的夾角

ρ——質量密度，kg/m3

σ——原子間平衡間距，?

χ——標度系數，用來調整范德華作用力強弱下角標

i，j，k——氧原子編號