基于COMSOL軟件的綠洲鹽漬化土壤中多離子耦合運移模型構建

焦會青，盛 鈺，趙成義，李保國 ※

（1．中國農業(yè)大學資源與環(huán)境學院，農業(yè)部華北耕地保育重點實驗室，國土資源部農用地質量與監(jiān)控重點實驗室，北京100193；2．中國科學院新疆生態(tài)與地理研究所，烏魯木齊 830011）

0 引 言

水資源短缺和土壤鹽漬化是制約干旱區(qū)農業(yè)發(fā)展的重要因素[1-2]，鹽堿土對作物生長的危害，不僅來自于耕層土壤較高的鹽分濃度，也來自于鹽分離子的組成差異[3]。土壤中離子之間相互影響、相互作用，由于各類鹽分離子的化學特性存在差異，使其在土壤中的遷移能力以及淋洗難易程度不同，對作物的影響機制和程度也有很大差異。因此，探求土壤中各鹽分離子的遷移規(guī)律可以為鹽漬化土壤的綜合治理以及高效利用提供科學依據[4]。

基于過程的數值模型是研究土壤水和溶質運移的有力工具[5-7]，Robbins等[8]將化學沉淀溶解和離子交換平衡模型子程序與一維非飽和水流耦合，并用蒸滲儀[9]和田間試驗[10]數據對模型進行了驗證。Chen等[11-12]建立了一維穩(wěn)態(tài)飽和流條件下Ca2+、Na+、SO42-和Cl-的耦合運移模型，以動力學模型描述硫酸鈣的沉淀溶解，以Gapon方程描述陽離子的交換過程。?im?nek和Suarez[13]、Suarez和?im?nek[14]建立了UNSATCHEM，用于模擬非飽和多孔介質中離子的平衡或者動力學反應過程，模型考慮溫度及CO2的影響，并且已經成功用于模擬CaSO4對鈉質土的改良[15]、蒸滲儀中滴灌條件下鹽分離子的運移[16]等過程，近些年，其中主要的離子反應模塊已經嵌入HYDRUS，但是由于缺乏相應的試驗數據，而且涉及參數較多，該模塊應用相對較少[17-18]。Lekakis和Antonopoulos[19]在一維土壤水氮運移模型WANISIM中加入離子（Ca2+、Na+、Mg2+）運移模塊，使模型可用于鹽分管理工作，但是模型中僅考慮了陽離子的交換和吸附。

COMSOL是一款高度集成的數值仿真軟件，基于有限元理論，以高效的計算性能和杰出的多場耦合分析能力實現任意多物理場的數值仿真。在地球科學領域，COMSOL可以模擬地下水流動、土壤中溶質運移和氣流形成規(guī)律、地熱、地面沉降等問題[20-24]，Rosenbom等[25]采用COMSOL評估了小尺度的空間變異對土壤中殺蟲劑降解和淋洗過程的影響；Wang等[26]采用COMSOL模擬了一維土柱中水鹽運移以及鹽分結晶過程，以此研究鹽分結晶對路基的影響程度；Mao等[27]基于COMSOL構建了土壤與地上部的耦合模型，模擬了殺蟲劑在土壤中遷移以及揮發(fā)至大氣的過程。COMSOL中提供了非飽和多孔介質溶質運移基本模塊，沒有考慮土壤中各組分之間的化學反應，因此，將COMSOL應用于非飽和多孔介質鹽分離子耦合運移的研究較少。但是，COMSOL支持自定義偏微分方程（PDE）的求解，用戶可以在該模式下自由定義所需方程并實現方程之間的耦合，而且在模型的構建中，材料屬性、源匯項以及各種邊界條件都可以使用任意類型的函數，比如分段、插值、解析等，這些函數既可以時間相關，也可以空間相關，其自變量既可以是獨立變量，也可以是求解變量本身[28]，借助于這些優(yōu)勢可以在COMSOL中進行物理問題的二次開發(fā)，而不需要另外編程求解，對于實現鹽分離子之間特定的化學反應有很大優(yōu)勢。

本研究采用COMSOL多孔介質和地下水流模塊模擬土壤水分運動，以溶質運移及鹽分離子反應理論為基礎，基于 COMSOL自定義偏微分方程模塊構建鹽漬化土壤SO42-、Ca2+、Na+、Cl-、Mg2+耦合運移模型，并以綠洲棉田膜下滴灌條件下主要鹽分離子的遷移為例對模型進行了驗證，為鹽分離子運移模擬研究提供了一種簡便的模型構建方法。

1 非飽和多孔介質鹽分多離子耦合運移數學模型

1.1 控制方程

土壤中水分運動以 Richards方程描述，其中，根系吸水基于潛在蒸騰量、根系分布以及土壤基質勢計算，具體方法見文獻[29]。

本研究主要探討 SO42-、Ca2+、Na+、Cl-、Mg2+的互相作用和遷移過程，離子的運移過程基于對流彌散方程構建，SO42-、Ca2+、Na+、Cl-、Mg2+的運移方程如下，

式中：θ為土壤體積含水量，L3/L3；c1-c5分別為土壤溶液中 SO42-、Ca2+、Na+、Cl-、Mg2+的濃度，N/L3；s2、s3、s5分別為Ca2+、Na+、Mg2+的吸附量，N/M；t為時間，T；xi為空間坐標，L；ui為土壤孔隙流速，L/T；ρb為土壤容重，M/L3；X為CaSO4的沉淀量，N/L3；Dij為水動力彌散系數，L2/T，計算方法參考文獻[30]。

1.2 離子反應

對于鹽漬化土壤，土壤膠體上供吸附陽離子的剩余電荷一般較少，因此吸附作用通常較小，主要發(fā)生的是陽離子的交換過程。如果土壤溶液中離子組成或濃度發(fā)生改變，勢必會影響各離子在溶液和膠體之間的平衡狀態(tài)。假設土壤中的陽離子交換量不變且等于交換性Ca2+、Mg2+、Na+之和，用 Gapon方程描述Ca2+、Na+、Mg2+之間的交換過程[29]，

式中：KNa/Ca為Na+和Ca2+的交換常數，(N/L3)-1/2；KMg/Ca為 Mg2+和 Ca2+的交換常數；γ2、γ3、γ5分別為 Ca2+、Na+、Mg2+的活度系數。

硫酸鈣的沉淀溶解滿足溶度積原理，反應以二級動力學過程描述[31]，

式中：Ksp為溶度積常數，(N/L3)2；kd為溶解常數，L3/(N·T)；kp為沉淀常數，L3/(N·T)；α(X)為單位躍遷函數，即 X ＞ 0時為1，X ≤ 0時為0。

1.3 活度系數

土壤溶液中可溶鹽濃度增加或減少會改變各離子的活度系數，進而影響陽離子交換和溶解沉淀過程，活度系數的估算可采用Davies方程[31]，但是當離子強度I ＞ 0.1 mol/L時，活度系數與離子強度不再是通用的函數關系，而是取決于特定離子的交互作用[32]，因此，這里也可以采用適用于特定條件下的函數方程。

2 基于COMSOL鹽分多離子耦合運移模型實現方法

以COMSOL多孔介質和地下水流模塊模擬土壤水分運動，而對于COMSOL沒有定義的物理問題，可以通過自定義偏微分方程組（PDEs）實現二次開發(fā)，COMSOL中提供的系數型PDE如下，

式中：u 為求解變量；ea、da、c、α、γ、β、a、f均為自定義系數。

所有自定義系數可根據實際物理問題定義為不同類型的函數，這些函數可以是不同階次的導數，既可以時間相關，也可以空間相關，其自變量既可以是獨立變量，也可以是求解變量本身。對于各離子的控制方程，首先根據公式（9）的形式進行整理，然后在相應的編輯框中定義各系數，對于沒有對應項的部分則可通過定義系數f實現，而不再需要編程求解。就離子運移方程而言，可設定為對流彌散項，f則可作為源匯項，表征離子間的化學反應。各離子運移方程通過土壤含水量和達西流速與 Richards方程耦合，又通過離子之間的化學反應實現各離子方程之間的耦合。

相應的，該模塊也提供了多種開放的邊界條件，比如Dirichlet邊界和通量邊界條件，通量邊界定義如下可構成方程的時變項，

式中，g和q均為自定義系數。

3 模型檢驗：綠洲棉田膜下滴灌土壤鹽分離子運移

以新疆綠洲膜下滴灌條件下的離子遷移過程為例，通過計算模擬值與實測值之間的平均絕對誤差（MAE）、均方根誤差（RMSE）、平均相對誤差（MRE）以及決定系數（R2）對模型進行檢驗。

3.1 試驗及模型設置

試驗在新疆阿克蘇水平衡試驗站開展，共設置 3個處理，單次灌水量分別為46.8 mm（T1）、39.6 mm（T2）、28.8 mm (T3)，在無降雨的情況下每隔6 d灌溉一次。棉花的種植模式如圖1所示，窄行10 cm，寬行65 cm，膜間距60 cm，一膜覆蓋四行棉花，滴灌管鋪設在寬行中間，滴頭間距30 cm。2015年和2016年膜下滴灌期間，在寬行和膜間的中間位置采集土壤樣品，采樣深度分別為0～20、20～40、40～60、60～80、80～100 cm，采用烘干法測定土壤含水量，每個灌水周期測定兩次。2015年8月1日至2015年8月25日在寬行和膜間深20、40、60、80 cm處埋設土壤溶液提取器，定點采集土壤溶液，共采集 6次。土壤溶液提取器主要由取樣器、負壓泵、采樣瓶組成，所用取樣器陶土頭的直徑為2 cm、長度為5 cm，當外加吸力大于土壤基質吸力時，使土壤溶液通過陶土頭進入采樣瓶，將收集的土壤溶液用于鹽分離子的測定。其中Cl-采用AgNO3滴定法測定，SO42-采用EDTA間接滴定法測定，Ca2+和Mg2+采用EDTA絡合滴定法測定，Na+和K+采用火焰光度計測定。灌溉水礦化度0.485 g/L，SO42-、Ca2+、Na+、Cl-、Mg2+濃度分別為 1.20、0.88、4.09、2.96、0.46 mmol/L。

模擬區(qū)域寬75 cm，高100 cm，模擬區(qū)域及水分運動的邊界條件見圖1，對于鹽分運動，考慮灌溉帶入的鹽分，但是由于降水量非常小，所以忽略降水帶入的鹽分，所有邊界均采用第三類邊界條件。將灌溉前的實測含水量和各離子含量進行二維插值，作為模擬初始值。

圖1 棉花種植模式及水分運動邊界條件示意圖Fig.1 s of planting mode of cotton and boundary conditions for soil water flow

土壤顆粒組成見表1，并根據Rosetta得到土壤水力學參數（van Genuchten方程[33]）的初始值，采用 2015年的土壤含水量數據對模型進行校準，并采用2016年的試驗數據驗證模型，率定的參數值見表1。計算彌散系數時所采用的縱向彌散度取12 cm，橫向彌散度為縱向彌散度的1/10，溶質在純水中的擴散系數取1.296 cm2/d，根系吸水的水分脅迫函數（Feddes模型）中h1、h2、h3、h4分別取-10、-25、-500、-14 000 cm[34]。

陽離子交換量CEC、Gapon常數KNa/Ca和KMg/Ca參考文獻[17，35]分別設定為 6.996 × 10-2mol/kg、0.50(mol/L)-1/2、0.52；硫酸鈣溶度積常數Ksp取2.399 × 10-5(mol/L)2[31]，硫酸鈣溶解速率與土壤水流速有關[36], kd取3.13 × 10-3L/(mol·s)，沉淀速率常數 kp取 2.22 × 10-3L/(mol·s)[35]。由于試驗地抽提的土壤溶液樣品中 90%的離子強度都大于0.1 mol/L，且土壤水溶性SO42?占全鹽比例約為 67%，Ca2+占全鹽比例約為 21%，假設土壤溶液中Ca2+和SO42?處于平衡狀態(tài)，根據CaSO4的溶度積常數和土壤溶液中 SO42?和 Ca2+的濃度反推出 SO42?和 Ca2+的活度系數，以離子強度為橫坐標，活度系數為縱坐標，繪制各個散點，并擬合得到指數方程，與Davies方程的對比結果見圖2。

表1 土壤顆粒組成及水力學參數Table 1 Particle size fractions and corresponding soil hydraulic parameter values

圖2 活度系數與離子強度的關系Fig.2 Relationship between activity coefficient and ionic strength

3.2 模擬結果驗證

表2給出了2015年和2016年各處理土壤含水量模擬值與實測值的統計結果，MAE值介于 0.023～0.033 cm3/cm3，RMSE 值介于0.030～0.040 cm3/cm3，與Phogat等[37]和Ramos等[18]得到的結果相似。

表2 土壤含水量模擬值與實測值的統計對比Table 2 Statistical comparison of simulated and measured soil water contents

2015年各離子模擬值與實測值統計對比結果見表3，Ca2+和SO42-的模擬效果總體上優(yōu)于Na+、Cl-和Mg2+，處理一和處理二優(yōu)于處理三，這是因為處理三灌溉量較小，只有少量的灌溉水能夠到達膜間位置，導致膜間始終保持在相對干燥的狀態(tài)，在整個觀測期內膜間深度 20 cm處僅收集到一次土壤溶液，而且為了收集到能夠滿足測量的水量，抽提時間相對較長，因此在膜間沒有合理的初始值導致了該處理模擬值與實測值較大的差異。通常情況下，RMSE ≥ MAE，而RMSE大于MAE的程度往往可以反映實測值與模擬值之間差異的變化程度，或者說在一定意義上可以反映最大差異的程度[18]，正因為處理三膜間位置較大的模擬誤差，導致處理三RMSE明顯大于MAE。該模型模擬值與實測值之間的差異，一方面可以歸因于土壤溶液離子濃度的測定方法，土壤溶液抽提器采集的是陶土頭周圍一定時間和空間范圍內的土壤溶液，該范圍與陶土頭周圍土壤的水力學性質和基質勢有關[38]，但是模型給出的是陶土頭相應位置點的數據；另一方面可歸因于模型參數的設置，陽離子交換和硫酸鈣沉淀溶解反應的相關參數均根據文獻給出，而這些參數與特定的試驗環(huán)境有關，也必然會導致模擬值與實測值的差異。獨立擾動各反應參數，當陽離子交換量 CEC、Na+和 Ca2+的交換常數 KNa/Ca、Mg2+和 Ca2+的交換常數KMg/Ca的變化幅度為±20%時，與本研究模擬結果相比，所有離子模擬值與實測值之間的 RMSE值分別波動-2.60%～3.99%、0.25%～3.45%、-1.45%～2.74%；當硫酸鈣沉淀常數 kp、硫酸鈣溶解常數 kd的變化幅度為±20%時，所有離子模擬值與實測值之間的RMSE值分別波動-0.29%～1.51%、-9.36%～11.10%，其中，kd對SO42-和Ca2+的模擬精度影響較大，當 kd的變化幅度為±10%時，其RMSE值波動范圍介于-4.69%～7.36%，這是因為灌溉使膜下土壤溶液中SO42-和Ca2+濃度明顯降低，打破了硫酸鈣的沉淀溶解平衡，主要發(fā)生著硫酸鈣的溶解過程，kd影響硫酸鈣的溶解速度，因此會對SO42-和Ca2+的模擬產生較大的影響。另外，氣溫對上部土層水溶性離子含量的變化有明顯的影響[39]，本研究沒有考慮溫度因素也可能是產生模擬誤差的原因之一。

3.3 膜下滴灌條件下各離子在土壤剖面上的運移特征

表3 土壤溶液離子濃度模擬值與實測值的統計對比Table 3 Statistical comparison of simulated and observed ion concentrations in soil solution

以處理二為例，圖3給出了一個灌水周期（7 d）內土壤含水量的變化過程，從第二天開始灌溉，共 24 h。滴灌結束后，滴灌管下方土壤含水量接近飽和，但是由于膜下滴灌條件下，僅有少量的灌溉水可以到達膜間位置，使膜間地表在灌溉后也保持在相對干燥的狀態(tài)，隨著時間的推移，滴灌頭周圍的土壤水分逐漸向深層和膜間運動，使膜下土壤含水量逐漸降低。由于土壤顆粒組成的差異以及地下水的補給等原因，32 cm以下土層的含水量總體上高于上層。

圖3 單個灌水周期內土壤含水量（cm3·cm-3）在剖面上的變化過程Fig.3 Change of soil water contents in profile during an irrigation cycle

圖 4給出了各離子在一個灌水周期內的變化過程，灌水結束后所有鹽分離子在膜下位置都明顯降低，主要是 0～40 cm土層，隨著時間的推移，該位置的離子濃度均有不同程度的提升。灌溉導致土壤溶液中Ca2+和SO42-濃度降低，打破硫酸鈣的沉淀溶解平衡，促使硫酸鈣溶解，直至土壤溶液中Ca2+和SO42-活度的乘積等于溶度積常數；Ca2+濃度的逐漸增加提高了其交換能力，因而進一步促進膠體上吸附的 Na+和 Mg2+釋放至土壤溶液；除此之外，根系的吸水過程也會導致離子濃度的升高，這也是 Cl-濃度升高的主要原因。Cl-、Na+在灌溉后濃度升高的幅度小于Ca2+和SO42-，淋洗效果較好，這與文獻[40,41]研究結果相似，由此看來，在該試驗條件下，滴灌不僅能降低根區(qū)土壤鹽分濃度，也能明顯降低對作物危害較大的Cl-和Na+的組成比例，為棉花的正常生長提供適宜的土壤環(huán)境。

膜間地表由于蒸發(fā)作用的影響，使各離子均在地表逐漸積累，從圖4可以看到，Cl-在膜間地表積累量最大，這主要是由于Cl-不參與化學反應，受土壤固相影響相對較小，移動性較強，而Ca2+和SO42-的活度乘積在大于溶度積常數時就會產生沉淀，使其維持在相對平衡的狀態(tài)。

3.4 活度系數對離子運移的影響

本研究采用指數方程（由試驗數據擬合得到）計算活度系數，為了分析活度系數對離子運移的影響，與基于Davies活度系數方程的模擬結果進行了對比，圖5給出了一個灌水周期內2種方法計算得到的活度系數分布。這兩種方法都是以離子強度為自變量，且隨著離子強度的增加而減小，因此整個剖面上活度系數的變化規(guī)律與離子濃度相反，即灌溉后膜下活度系數增加，隨著時間的推移又逐漸減小，而在膜間近地表隨著蒸發(fā)積鹽的進行有小幅度的減小。同時，也可以看到在整個土壤剖面上，指數方程計算得到的活度系數小于Davies方程，例如在灌溉后，在整個剖面上，指數方程計算得到的活度系數介于0.07～0.37，而Davies方程計算得到的活度系數介于 0.29～0.57。

圖4 單個灌水周期內土壤溶液離子濃度在剖面上的變化過程Fig.4 Changes of various ion concentrations in soil solution in profile during an irrigation cycle

圖6為基于Davies活度系數估算方法得到的各離子的模擬結果，與基于指數方程得到的模擬結果（圖4）相比，Ca2+和 SO42-的差異最明顯。在第一天尚未灌溉時，整個剖面上土壤溶液中Ca2+和SO42-濃度小于圖4中模擬結果，這是因為Davies方程計算的活度系數大于指數方程，使土壤溶液中Ca2+和SO42-活度乘積大于溶度積常數，促進了 CaSO4的沉淀過程。在第二天灌溉后，膜下淺層土壤中Ca2+和SO42-都降低，此時CaSO4開始溶解，其溶解速度與溶度積常數和這兩種離子活度乘積的差異成正比，因此活度系數較高時 CaSO4的溶解速度較低，導致在灌溉以后膜下Ca2+和SO42-的濃度依然低于圖4中的模擬結果。

圖5 單個灌水周期內活度系數在剖面上的變化過程Fig.5 Changes of activity coefficient in the profile during an irrigation cycle

兩種活度系數計算方法在膜間地表產生的差異則更為明顯，基于指數方程得到的Ca2+和SO42-濃度在逐漸升高，而基于Davies方程得到的Ca2+和SO42-濃度有小幅度的減小，這是因為Davies方程計算得到的活度系數較大，模擬開始時刻該處Ca2+和SO42-的活度乘積大于溶度積常數，所以在整個灌水周期內都處于沉淀過程，使 Ca2+和SO42-濃度逐漸降低；而指數方程計算的活度系數相對較小，兩種離子活度的乘積尚未達到溶度積常數，因此使離子在此聚積而尚未沉淀。

Na+和Mg2+的變化趨勢與圖4結果類似，但是總體上略小于圖4，這是因為Na+、Mg2+和Ca2+處于動態(tài)的交換過程中，當土壤溶液中 Ca2+濃度較小時，其交換能力較小，土壤溶液中的Na+和Mg2+則更易被土壤膠體吸附。基于Davies活度系數估算方法得到的Ca2+、SO42-、Na+、Mg2+的模擬結果與實測值間的 RMSE值分別為 7.87、16.81、6.55、4.93 mmol/L，明顯大于基于指數方程得到的結果。因此，活度系數的準確估算對模擬結果的準確性有重要影響，尤其是參與沉淀溶解反應的Ca2+和SO42-，當離子強度I ＞ 0.1mol/L時，活度系數與離子強度不再是通用的函數關系[32]，采用通用的函數關系可能會帶來較大的模擬誤差。

圖6 基于Davies方程單個灌水周期內土壤溶液離子濃度在剖面上的變化Fig. 6 Changes of ion concentrations of soil solution in profile during an irrigation cycle based on Davies Equation

4 結論

1）COMSOL在材料屬性、源匯項以及邊界條件的設定上非常靈活，而且支持自定義偏微分方程的求解，基于這些優(yōu)勢，在多孔介質和地下水流模塊模擬非飽和土壤水流的基礎上，自定義偏微分方程組構建鹽漬化土壤SO42-、Ca2+、Na+、Cl-、Mg2+耦合運移模型，考慮陽離子交換過程以及硫酸鈣的沉淀溶解反應，并以綠洲棉田膜下滴灌田間試驗對模型進行了驗證，結果表明該模型能夠較好地反映土壤中鹽分離子的動態(tài)變化規(guī)律。但是，該模型未考慮碳酸鈣的沉淀溶解過程，因此，對于碳酸鹽質量分數較高的土壤，該模型不適用，可以進一步在COMSOL自定義偏微分方程模塊下定義該過程。

2）膜下滴灌條件下，膜下土壤中的鹽分離子有不同程度的淋洗，其中，Cl-、Na+在灌溉后濃度升高速度小于Ca2+和 SO42-，淋洗效果較好。膜間地表由于蒸發(fā)作用的影響，使各離子在地表逐漸積累，Cl-移動性較強，在膜間地表積累程度最大。

3）活度系數估算方法對模擬結果的準確性有重要影響，尤其是鹽分含量較高時，采用通用的函數關系可能會帶來較大的模擬誤差。