(1.長(zhǎng)江水利委員會(huì)水文局 長(zhǎng)江下游水文水資源勘測(cè)局，江蘇 南京 210011； 2. 南京水利科學(xué)研究院 水文水資源與水利工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210029)

我國(guó)巖溶分布面積廣，引起的資源環(huán)境及工程問(wèn)題也較為突出，如巖溶區(qū)隧洞涌突水、水庫(kù)滲漏、邊坡穩(wěn)定等問(wèn)題，而此類問(wèn)題多與巖溶演化發(fā)育有關(guān)[1-4]。地下水在運(yùn)移過(guò)程中與可溶性碳酸鹽類發(fā)生反應(yīng)，引起礦物溶解，在以裂隙為主的巖溶系統(tǒng)中使得裂隙開度發(fā)生擴(kuò)張，并導(dǎo)致滲流場(chǎng)、化學(xué)場(chǎng)等發(fā)生不可逆的改變，此過(guò)程可通過(guò)滲流-溶解耦合模型加以反映刻畫[5]。模型主要由水流模塊、溶質(zhì)運(yùn)移模塊、裂隙開度時(shí)變模塊等組成，并通過(guò)各模塊變量間的函數(shù)關(guān)系實(shí)現(xiàn)模型耦合[6-7]。其中，礦物溶蝕速率作為模型的重要組成，在以上3個(gè)模塊中均有所體現(xiàn)。已有研究表明，在不考慮CO2對(duì)反應(yīng)速率影響的封閉環(huán)境中，溶蝕速率主要受表面反應(yīng)和擴(kuò)散遷移雙重因素的聯(lián)合控制，該兩項(xiàng)進(jìn)程逐次進(jìn)行，溶蝕速率值為兩者間相對(duì)較小者[8]。但在實(shí)際模擬應(yīng)用中，由于該速率方程表達(dá)較為復(fù)雜，常僅考慮受單因素控制，即只受表面反應(yīng)或擴(kuò)散遷移其中某一因素控制，從而使得不同研究中溶蝕速率方程有所差異[9-11]?？梢钥闯?，以往研究多側(cè)重于研究裂隙系統(tǒng)的發(fā)育演變，以及由此帶來(lái)的滲流動(dòng)態(tài)、溶液組分等的改變，而未涉及在模擬過(guò)程中，當(dāng)溶蝕速率表達(dá)不同時(shí)，對(duì)模擬結(jié)果存在何等影響。

針對(duì)目前研究中存在的不足，本文旨在量化解析不同溶蝕速率方程對(duì)模型結(jié)果的影響。首先根據(jù)溶蝕方程機(jī)理的不同，分為單因素控制(表面反應(yīng)或擴(kuò)散遷移)以及兩者聯(lián)合控制。然后假設(shè)裂隙形態(tài)皆為平行板裂隙，分別模擬獲得不同溶蝕速率下裂隙開度、滲流性態(tài)以及溶液組分等的空間分布和時(shí)間演變規(guī)律，并進(jìn)行比較獲得其異同。研究成果可為今后實(shí)際應(yīng)用中根據(jù)需要選擇更為合適的模型提供指導(dǎo)。

1 溶蝕機(jī)理與速率模型

灰?guī)r單裂隙滲流-溶解耦合模型主要由水流模塊、溶質(zhì)運(yùn)移模塊和裂隙開度時(shí)變模塊3部分組成，具體可見(jiàn)文獻(xiàn)[7，10]，在此不再贅述。由于灰?guī)r溶蝕速率是滲流-溶解耦合模型的核心內(nèi)容之一，在模型3個(gè)模塊中皆有涉及。根據(jù)其受控機(jī)理的不同，一般包括以下3種：① 僅受表面反應(yīng)單一因素控制；② 僅受擴(kuò)散遷移單一因素控制；③ 受表面反應(yīng)和擴(kuò)散遷移雙重因素聯(lián)合控制。以上各機(jī)理及速率模型如下所述。

(1) 表面反應(yīng)控制，即僅受可溶性礦物的溶解動(dòng)力學(xué)控制，其速率大小與溶液中對(duì)應(yīng)礦物的相對(duì)飽和度有關(guān)。對(duì)于碳酸鹽類，在接近平衡態(tài)和遠(yuǎn)離平衡態(tài)時(shí)，其速率表達(dá)略有不同。當(dāng)遠(yuǎn)離平衡態(tài)時(shí)反應(yīng)相對(duì)較快，而接近飽和平衡態(tài)時(shí)，其速率陡降，該過(guò)程可表示為

(1)

式中，k1、k2為動(dòng)力學(xué)反應(yīng)常數(shù)；n1、n2為指數(shù)，通常n1為等于1或?yàn)榻咏?的數(shù)，而n2在3～6之間；Cs為拐點(diǎn)濃度值，一般為0.7Ceq～0.9Ceq之間。

(2) 擴(kuò)散遷移控制，即對(duì)應(yīng)離子組分從固相表面進(jìn)入溶液中的遷移速率，一般認(rèn)為該速率不僅與溶液中離子濃度有關(guān)，還與該位置處裂隙開度b有關(guān)，表達(dá)如下：

(2)

式中，Dm為離子在水溶液中的擴(kuò)散系數(shù)；ShD為Sherwood數(shù)，在層流情形下為8.24[12]；Ceq為離子飽和濃度。

需要指出的是，國(guó)內(nèi)在模擬滲流-溶解耦合問(wèn)題時(shí)主要基于“薄膜假設(shè)”，認(rèn)為緊貼裂隙壁有一層極薄的“薄膜”，內(nèi)部溶液濃度飽和，薄膜內(nèi)組分可進(jìn)入擴(kuò)散區(qū)，導(dǎo)致薄膜內(nèi)溶液濃度降低；此時(shí)，固壁表面發(fā)生溶解作用，離子進(jìn)入“薄膜”補(bǔ)充其內(nèi)物質(zhì)損失，使之重新達(dá)到飽和；同時(shí)，由于固壁表面物質(zhì)的溶失，造成裂隙開度的增大[13]。該假設(shè)本質(zhì)上也為擴(kuò)散遷移控制，但在該速率模型中不考慮裂隙開度變化對(duì)溶蝕速率的影響。

(3) 表面反應(yīng)和擴(kuò)散遷移兩因素聯(lián)合控制，即認(rèn)為在灰?guī)r溶解過(guò)程中，表面反應(yīng)與擴(kuò)散過(guò)程逐次依序進(jìn)行，因此總的溶解速率取決于各時(shí)刻以上兩過(guò)程中的最慢者，即以上式(1)與式(2)計(jì)算結(jié)果間相對(duì)較小者作為該時(shí)刻下灰?guī)r的溶解速率。

2 不同溶蝕速率對(duì)模型影響分析

2.1 模型簡(jiǎn)介

巖溶地區(qū)裂隙相互交錯(cuò)形成裂隙網(wǎng)絡(luò)，其真實(shí)形態(tài)極為復(fù)雜，因此常以對(duì)單裂隙滲流特性的研究作為認(rèn)識(shí)巖體裂隙整體滲流特征的基礎(chǔ)。將裂隙壁簡(jiǎn)化為光滑平板，同時(shí)隙寬相對(duì)較小時(shí)，可采用立方定律研究模型隙寬與單寬流量之間的關(guān)系，在此基礎(chǔ)上加以修正可推廣至粗糙裂隙及裂隙網(wǎng)絡(luò)。因此，平行板裂隙下滲流特性的研究作為真實(shí)裂隙網(wǎng)絡(luò)的基礎(chǔ)具有較為重要的意義[14-15]。

考慮研究區(qū)為一長(zhǎng)為L(zhǎng)、寬為W的單裂隙，其初始裂隙隙寬為b0，裂隙兩端水頭分別為Hin和Hout，流入溶液的濃度為Cin，裂隙內(nèi)初始水頭和初始濃度分別為H0和C0，模型示意圖如圖1所示，具體取值見(jiàn)表1。

圖1 模型示意Fig.1 Sketch of model

L/mW/mαLαTM/(g·mol-1)ρ/(kg·m-3)Cin/(mol·m-3)C0/(mol·m-3)b0/mmH0/mHin/mHout/m10.50.050.0210026600.0500.200.10

2.2 溶蝕速率參數(shù)及比較

表面反應(yīng)和擴(kuò)散遷移控制中涉及到的主要參數(shù)如表2所示。為研究不同溶蝕速率方程對(duì)模型結(jié)果的影響，分別建立以下模型，其溶蝕速率對(duì)應(yīng)控制因素為：模型1為僅表面反應(yīng)單一因素控制，模型2為僅擴(kuò)散遷移單一因素控制，模型3為表面反應(yīng)及擴(kuò)散遷移兩者聯(lián)合控制。

表2 反應(yīng)速率方程中相關(guān)參數(shù)Tab.2 Parameters in reaction rate equation

注：表中參數(shù)引自文獻(xiàn)[8]。

在有限元軟件COMSOL Multiphysics軟件中建立尺寸為L(zhǎng)×W(長(zhǎng)×寬)的二維研究區(qū)域，并設(shè)定相應(yīng)的模型方程和定解條件。采用瞬態(tài)分離式求解器( time dependent segregated) 對(duì)以上3個(gè)模型分別進(jìn)行求解。在其它條件相同情形下，模擬時(shí)段(2 000 d)內(nèi)選取500,1 000,1 500 d和2 000 d，在上述時(shí)刻時(shí)各模型裂隙開度變化如圖2所示。由圖2可以看出：

(1) 3個(gè)模型在上游進(jìn)水口位置溶蝕皆較為嚴(yán)重，而向下游溶蝕逐漸減緩，從溶蝕鋒面來(lái)看，在前期(500 d和1 000 d)皆較為平直，而在中后期(1 500 d和2 000 d)出現(xiàn)溶蝕差異不均衡現(xiàn)象。

(2) 3個(gè)模型的不同主要體現(xiàn)在溶蝕“程度”和“深度”兩方面。從溶蝕“程度”來(lái)看，在研究區(qū)上游側(cè)，尤其是靠近進(jìn)水口部位，模型1結(jié)果較模型2和模型3明顯偏大，而在研究區(qū)的中后部，模型1和模型2未發(fā)生明顯溶蝕，其開度b較模型3明顯偏??；從溶蝕“深度”看，模型3>模型2>模型1，其中模型3于2 000 d時(shí)在研究區(qū)靠近軸部的部位形成了一條較寬的集中滲流通道，且基本貫穿整條裂隙，模型2則形成了局部的滲流通道，達(dá)到了裂隙中部，但未貫穿，而模型1的溶蝕僅發(fā)生在靠近進(jìn)水口的較短距離內(nèi)(<0.2 m)。

圖2 不同時(shí)刻下裂隙開度b的分布(由上至下依次為500,1 000 ,1 500 d和2 000 d)Fig.2 Distribution of fracture aperture at different time(from top to bottom are in 500,1 000,1 500 d and 2 000 d)

為更直觀地比較3個(gè)模型各時(shí)刻下開度在沿裂隙延伸方向上的變化，以y=0.25 m作為典型斷面，以上時(shí)刻各模型開度b在該典型斷面上的分布見(jiàn)圖3。

由圖3可以看出：相同時(shí)刻下，在裂隙前端靠近進(jìn)水口部位(<0.1 m)，模型2和模型3的開度分布較為近似，且都明顯小于相同位置處模型1的。這表明模擬時(shí)段較長(zhǎng)時(shí)，該部位溶解主要受擴(kuò)散遷移控制，而采用表面反應(yīng)控制的動(dòng)力學(xué)表達(dá)，會(huì)導(dǎo)致該部位灰?guī)r溶蝕的過(guò)量估計(jì)。在研究區(qū)中后段(>0.4 m)，尤其是模擬時(shí)刻的后期(如2 000 d)，模型1和模型2的結(jié)果明顯小于模型3的，表明單一因素控制的模型在研究區(qū)范圍相對(duì)較大時(shí)，其位于研究區(qū)中后部的結(jié)果誤差相對(duì)較大。

可通過(guò)計(jì)算獲得模擬時(shí)段內(nèi)研究區(qū)平均隙寬隨時(shí)間的變化。平均隙寬一方面可反映裂隙形態(tài)的變化，另一方面也可反映灰?guī)r裂隙整個(gè)溶失量的變化。同時(shí)，為更直觀地比較各模型間的差異，以模型3的結(jié)果作為基準(zhǔn)，計(jì)算獲得模型1、模型2與模型3間的比例。各模型平均隙寬及對(duì)應(yīng)比例隨時(shí)間變化見(jiàn)圖4。

由圖4可以看出，3個(gè)模型計(jì)算獲得的實(shí)際平均隙寬在整個(gè)模擬時(shí)段內(nèi)皆呈增大的趨勢(shì)，但彼此間的大小關(guān)系在不同時(shí)刻有所不同。在該算例中，各模型間的差異性變化可分為以下4個(gè)階段：① 300 d之前，各模型的平均隙寬表現(xiàn)為模型2>模型1>模型3；② 300～510 d之間，表現(xiàn)為模型1>模型2>模型3；③ 在510～1 370 d之間，則為模型1>模型3>模型2；④ 在1 370 d之后，模型3>模型1>模型2，且這段期間模型3與其他兩個(gè)模型間的差異越來(lái)越大，在模擬結(jié)束時(shí)模型1和模型2僅為模型3值的0.739和0.593。以上4個(gè)階段可簡(jiǎn)單概括為：前期單一因素控制下模型平均隙寬變化大于兩因素聯(lián)合控制的；而后模型3(兩因素聯(lián)合控制)的隙寬變化依次超過(guò)模型2(擴(kuò)散遷移控制)和模型1(表面反應(yīng)控制)，且與后兩者間的差異越來(lái)越顯著。

由以上可知，不同溶蝕速率方程對(duì)裂隙開度的變化有所影響，因此通過(guò)整個(gè)裂隙的流量也會(huì)有所差異。同樣以模型3的結(jié)果作為基準(zhǔn)，分別計(jì)算獲得模型1、模型2與模型3間的比例，則出口位置處各模型通過(guò)裂隙的流量Q及對(duì)應(yīng)比例隨時(shí)間的變化如圖5所示。

由圖5可以看出：① 3種模型中通過(guò)裂隙的流量，總體上隨時(shí)間的推移皆呈逐漸增長(zhǎng)的趨勢(shì)，而相同時(shí)刻下，三者比較，總體上呈“模型3>模型1>模型2”。該結(jié)果表明，由單一因素控制溶解速率的模型計(jì)算獲得的流量結(jié)果普遍低于實(shí)際通過(guò)裂隙的流量，其中尤以模型2的差別相對(duì)較大。在290d時(shí)，模型2與模型3的差別即超過(guò)了10%，而模型1在1080d后才達(dá)到該數(shù)值。②模型3中流量增長(zhǎng)出現(xiàn)“臨界點(diǎn)”，在該時(shí)刻前增長(zhǎng)速率較為平緩，而該時(shí)刻后，其速率陡增，并最終達(dá)到初始時(shí)的96.57倍，此與模型3中裂隙形成貫穿型滲流通道有關(guān)；而模型1和模型2在整個(gè)模擬時(shí)段內(nèi)流量皆呈緩慢增長(zhǎng)的趨勢(shì)，無(wú)明顯的流量激增時(shí)刻，而終了時(shí)刻時(shí)的流量也僅為初始時(shí)的3.22倍和1.83倍。以上結(jié)果表明，在對(duì)通過(guò)裂隙的滲流量進(jìn)行預(yù)測(cè)時(shí)，采用單因素控制的模型會(huì)低估通過(guò)裂隙的滲流量，尤以擴(kuò)散遷移控制的模型表現(xiàn)明顯。

圖3 不同時(shí)刻下y=0.25 m處裂隙開度變化Fig.3 Fracture aperture along y=0.25 m

圖6 不同時(shí)刻下研究區(qū)內(nèi)Ca2+濃度的分布(由上至下依次為500，1 000，1 500 d和2 000 d)Fig.6 Distribution of Ca2+ concentration in 500,1 000,1 500 d and 2 000 d

圖4 不同模型裂隙平均隙寬及相應(yīng)比例Fig.4 Average aperture width and corresponding proportion of different models

圖5 不同模型通過(guò)裂隙流量及比例過(guò)程線Fig.5 Flow through aperture of different models and corresponding proportion curve

除裂隙開度、滲流場(chǎng)外，模型間化學(xué)場(chǎng)分布亦有所不同，其中Ca2+分布如圖6所示。由圖6可以看出各模型計(jì)算獲得的Ca2+分布具有以下異同點(diǎn)。

(1)對(duì)于每個(gè)模型而言，Ca2+濃度沿裂隙延伸方向逐漸增大，在進(jìn)水口位置處等于流入溶液的Ca2+濃度，而在一段距離后接近飽和；且與裂隙開度分布相似，鋒面在前期較為平直，后在中后期出現(xiàn)溶蝕差異。

圖7 不同時(shí)刻y=0.25 m處Ca2+濃度分布Fig.7 Concentration of Ca2+ along y=0.25 m

(2)模型3由于在模擬時(shí)段后期在近裂隙軸部部位形成較為集中的滲流通道，且貫穿整個(gè)裂隙，因此在該段時(shí)期裂隙溶液中Ca2+濃度相對(duì)較小，而相同時(shí)段內(nèi)模型1和模型2該部位的Ca2+濃度仍相對(duì)較大。

(3)在模擬時(shí)段前期(<1 500 d)，盡管3個(gè)模型在裂隙中后部的溶液Ca2+濃度皆接近飽和，但其值大小仍有區(qū)別，其中模型2的濃度值已基本等于灰?guī)r溶解的飽和濃度值，即0.1768 mol/m3；而模型1和模型3的最大值仍小于0.16 mol/m3，這與表面反應(yīng)控制在溶液達(dá)到拐點(diǎn)濃度后，其反應(yīng)速率大大減慢有關(guān)。

為更清晰地比較各模型中Ca2+濃度沿裂隙延伸方向上的變化，以y=0.25 m作為典型斷面，則各模型在該斷面上Ca2+濃度分布如圖7所示。由圖7可以看出：

(1)無(wú)論是單一因素控制的模型1和模型2，還是雙重因素控制的模型3，Ca2+濃度在沿裂隙延伸方向上皆呈逐漸增大的趨勢(shì)，且在裂隙前部增長(zhǎng)較快，而在裂隙中后部增長(zhǎng)緩慢。

(2)模型1與模型3在裂隙中后部Ca2+的濃度值呈緩慢增長(zhǎng)的趨勢(shì)，但仍未達(dá)到飽和狀態(tài)，而模型2則達(dá)到飽和平衡態(tài)。

3 結(jié) 論

在建立灰?guī)r單裂隙滲流-溶解耦合模型的基礎(chǔ)上，研究了不同溶蝕速率對(duì)模型結(jié)果的影響，并以一平行板裂隙中滲流溶解為例，得出如下結(jié)論。

(1) 對(duì)于裂隙開度而言，從溶蝕深度來(lái)看，模擬前期3種模型基本一致，而中后期雙重因素下模型明顯更深，從而形成貫穿裂隙的滲流通道；就溶蝕程度來(lái)看，3種模型下進(jìn)水口位置溶蝕更為嚴(yán)重，其中尤以表面反應(yīng)單因素控制的更為明顯，其在上游位置的溶蝕明顯大于另外兩種模型。

(2) 就總的溶失量而言，前期單一因素控制下模型平均隙寬變化大于兩因素聯(lián)合控制的，但相差不大；而后兩因素聯(lián)合控制依次超過(guò)擴(kuò)散遷移控制和表面反應(yīng)控制，且與后兩者間的差異越來(lái)越顯著。

(3) 從通過(guò)裂隙的滲流量來(lái)看，采用單因素控制的模型會(huì)低估通過(guò)裂隙的滲流量，尤以擴(kuò)散遷移控制的模型表現(xiàn)明顯。同時(shí)，溶蝕使得裂隙中可形成貫穿型集中滲流通道，從而引起流量陡增，單因素控制的模型預(yù)測(cè)的該臨界點(diǎn)遠(yuǎn)遠(yuǎn)晚于實(shí)際。

(4) 從化學(xué)場(chǎng)分布來(lái)看，模擬前期3種模型差別并不顯著，但隨著兩因素聯(lián)合控制的模型形成貫穿型滲流通道，其裂隙內(nèi)Ca2+濃度分布明顯小于另外另種模型。此外，對(duì)于擴(kuò)散遷移控制的模型，其在裂隙下游側(cè)溶液中Ca2+濃度達(dá)到飽和，其值明顯大于另兩種模型。

(5) 當(dāng)采用單因素控制對(duì)溶解速率進(jìn)行簡(jiǎn)化時(shí)，會(huì)對(duì)模型結(jié)果造成一定影響，但可根據(jù)應(yīng)用中關(guān)注的內(nèi)容(如裂隙開度、溶蝕量、滲流量、Ca2+濃度等)以及模擬時(shí)段的長(zhǎng)度，從表面反應(yīng)和擴(kuò)散遷移中選擇合適的模型，亦可滿足一定的精度要求。

需要指出的是，本文的研究主要基于灰?guī)r單裂隙，盡管單裂隙溶蝕擴(kuò)展特征及規(guī)律是巖溶發(fā)育演化研究的前提和基礎(chǔ)，但其與實(shí)際的巖溶發(fā)育演化有一定的差別。因此，研究不同溶蝕速率對(duì)粗糙裂隙形態(tài)以及裂隙網(wǎng)絡(luò)中滲流-溶解耦合模型的影響是下一步研究的重點(diǎn)。