，，，，，

(1.南京大學(xué)地球科學(xué)與工程學(xué)院水科學(xué)系/教育部表生地球化學(xué)重點實驗室，江蘇 南京 210023；2.南京師范大學(xué)虛擬地理環(huán)境教育部重點實驗室，江蘇 南京 210046)

巖溶管道中懸浮物運(yùn)移及等參替代數(shù)值模型研究
——以廣西某排泥庫懸浮物運(yùn)移預(yù)測為例

周宇超1，祝曉彬1，吳吉春1，吳劍鋒1，曾獻(xiàn)奎1，郭飛2

(1.南京大學(xué)地球科學(xué)與工程學(xué)院水科學(xué)系/教育部表生地球化學(xué)重點實驗室，江蘇 南京 210023；2.南京師范大學(xué)虛擬地理環(huán)境教育部重點實驗室，江蘇 南京 210046)

暗河發(fā)育巖溶地區(qū)地下暗河中污染物的運(yùn)移模擬研究一直是難點。實際工程應(yīng)用中，由于傳統(tǒng)數(shù)值模擬模型沒有參數(shù)反應(yīng)懸浮物的沉降作用，通常僅概化考慮懸浮物在地下暗河中的遷移而采用河流一維穩(wěn)態(tài)水質(zhì)模型進(jìn)行預(yù)測，缺少其它方法進(jìn)行有效驗證。本文通過實驗和數(shù)值模擬相結(jié)合的方法，研究了管徑、懸浮物初始濃度、粒徑、流速對管道中懸浮物沉降系數(shù)的影響，并初步探究數(shù)值模型機(jī)理相近的吸附系數(shù)等效替代一維穩(wěn)態(tài)水質(zhì)模型沉降系數(shù)的可行性，再進(jìn)一步根據(jù)廣西某巖溶區(qū)排泥庫工程構(gòu)建具體條件下吸附系數(shù)與沉降系數(shù)的等參替代關(guān)系，最后將其應(yīng)用于數(shù)值模型模擬。計算結(jié)果表明經(jīng)過參數(shù)等效替代后的數(shù)值模擬方法比一維穩(wěn)態(tài)水質(zhì)模型更為合理。

巖溶；一維穩(wěn)態(tài)水質(zhì)模型；數(shù)值模擬模型；吸附系數(shù)；沉降系數(shù)

我國巖溶區(qū)面積占國土面積近1/3，近年來，隨著巖溶區(qū)經(jīng)濟(jì)的持續(xù)發(fā)展和礦山的大量開采利用，過程中產(chǎn)生的洗礦泥漿等導(dǎo)致的地下水環(huán)境問題受到越來越廣泛的關(guān)注。有效保護(hù)巖溶區(qū)地下水環(huán)境的前提是對地下水流和溶質(zhì)運(yùn)移規(guī)律的掌握。

而巖溶區(qū)，特別是暗河發(fā)育地區(qū)，含水介質(zhì)呈高度管道化和強(qiáng)烈非均質(zhì)性，使得該區(qū)域地下水中污染物運(yùn)移模擬的研究一直十分困難。實際工程中常用的一維河流穩(wěn)態(tài)水質(zhì)模型僅將暗河管道作為研究對像，主要考慮污染物的初始濃度、流速和污染物的沉降系數(shù)等，而不考慮暗河管道之外區(qū)域?qū)ξ廴疚镞\(yùn)移的影響，其所得結(jié)果為一粗略近似解。對于不同分布式模型[1-4]，模型考慮的因素越多，要求獲取的巖溶含水系統(tǒng)特征和水文地質(zhì)參數(shù)越精細(xì)、越多。而實際巖溶含水系統(tǒng)的暗河發(fā)育位置、巖溶管道結(jié)構(gòu)、水流特征、管道-裂隙交互參數(shù)等的獲取均很難，涉及參數(shù)越多，模型實際應(yīng)用越困難。因此，實際工程應(yīng)用中，理論模型與實際條件之間達(dá)到一定平衡是關(guān)鍵。分布式模型中，等效多孔介質(zhì)模型[5]因其研究最為成熟、包含此模型的計算軟件也最多，涉及參數(shù)相對較少，操作簡便且可以較好處理巖溶含水介質(zhì)的非均勻性，被廣泛用于模擬巖溶暗河管道的溶質(zhì)運(yùn)移問題。

等效多孔介質(zhì)模型模擬預(yù)測巖溶地下暗河懸浮物運(yùn)移時，目前常用數(shù)值模擬軟件Feflow由于缺少反映懸浮物運(yùn)移過程中的沉降特征參數(shù)而使其應(yīng)用受到限制。本文通過實驗研究特定條件下數(shù)值模型吸附系數(shù)與河流解析模型沉降系數(shù)的等參替代關(guān)系，并以地下暗河發(fā)育的廣西某巖溶區(qū)排泥庫工程主要污染物——懸浮物SS在事故工況下的運(yùn)移預(yù)測為例，將等參代替數(shù)值模擬模型用于模擬巖溶地下暗河中懸浮物的運(yùn)移，為數(shù)值方法在巖溶區(qū)地下水環(huán)境影響預(yù)測評價中的應(yīng)用提供成功實踐。

1 巖溶暗河懸浮物運(yùn)移預(yù)測評價方法

1.1 河流一維穩(wěn)態(tài)水質(zhì)模型

河流水質(zhì)模型是巖溶管道發(fā)育地區(qū)溶質(zhì)運(yùn)移模擬最常用方法[6]。由于泥漿中污染物僅為泥質(zhì)懸浮物SS，作為保守性污染物，不考慮遷移過程中污染物的源匯項和污染物在地下暗河的吸附、揮發(fā)、生物化學(xué)反應(yīng)等，適用于河流一維穩(wěn)態(tài)水質(zhì)模型[7]：

C=C0×exp[-Kx/(86400u)]

(1)

式中：C0為初始濃度，mg/L；K為沉降系數(shù)，1/d；u為地下河流速，m/s；x為沿地下河方向距離，m；c為位于污染源(排放口)下游x處的水質(zhì)濃度，mg/L；

其中，初始濃度采用充分混合模式求?。?/p>

(2)

式中：C0為初始濃度，mg/L；Qp為污水排放量，m3/s；Cp為污染物排放濃度， mg/L；Qh為上游來水量，m3/s；Ch為上游來水污染物濃度，mg/L。

1.2 數(shù)值方法

根據(jù)研究區(qū)水文地質(zhì)條件，通過分析地下水補(bǔ)、徑、排特征，將研究區(qū)的地下水流概化成非均質(zhì)、各向異性、三維非穩(wěn)定潛水地下水流系統(tǒng)[7]，對應(yīng)的數(shù)學(xué)模型如下：

1.2.1 地下水流數(shù)學(xué)模型

(3)

式中：H0為水頭初值；φ1為第一類邊界Γ1上的已知函數(shù)；q為第二類邊界Γ2上的單位面積側(cè)向補(bǔ)給量。n為邊界Γ2的外法線方向；Ki,j為各方向上的滲透系數(shù)，ω為源匯項，包括降水入滲補(bǔ)給、蒸發(fā)、井的抽水量和泉的排泄量等。

基于研究區(qū)復(fù)雜邊界及地下暗河管道的刻畫，本文選用有限元法。利用FEFLOW求解水流模型，考慮到模型特點，選擇Galerkin有限元法求解，對于潛水面的處理，選擇迭代法求解。

1.2.2 溶質(zhì)運(yùn)移數(shù)學(xué)模型

懸浮物在地下水中的運(yùn)移為三維水流作用下的對流彌散問題，溶液密度假定不變，考慮局部平衡吸附和一級不可逆動力反應(yīng)。溶質(zhì)運(yùn)移的三維水動力彌散方程的數(shù)學(xué)模型如下[8]：

(4)

式中：

C為為地下水中組分的溶解相濃度，ML-3；θ為為地層介質(zhì)的孔隙度，無量綱；t為為時間，T；xi為沿直角坐標(biāo)系軸向的距離，L；Dij為水動力彌散系數(shù)張量，L2T-1；Vi為孔隙水平均實際流速，LT-1；qs為單位體積含水層流量，代表源匯項，L3T-1；Cs為源匯項中組分的濃度，ML-3；∑Rn為化學(xué)反應(yīng)項，ML-3T-1。

1)初始條件

將補(bǔ)給濃度邊界和注水井處的初始濃度定為C0，其余地方均為0 mg/L，具體表述為：

2)邊界條件

本次模擬將含水層各個邊界均看做二類邊界條件(Neumann邊界)，且穿越邊界的彌散通量為0，具體可表述為：

式∞=中：Г2為Neumann邊界

溶質(zhì)運(yùn)移模型采用FEFLOW求解。

2 巖溶管道懸浮物運(yùn)移等參替代研究

研究基于不同條件管道懸浮物運(yùn)移室內(nèi)實驗以及地下水等效介質(zhì)數(shù)值模擬模型和河流一維穩(wěn)態(tài)水質(zhì)模型。首先設(shè)置不同實驗條件，通過室內(nèi)實驗獲取對應(yīng)懸浮物濃度。然后以實驗設(shè)定條件為數(shù)值模擬模型和一維穩(wěn)態(tài)模型的構(gòu)建基礎(chǔ)，通過擬合實驗實測數(shù)據(jù)，分別求得數(shù)值模型中的吸附系數(shù)和一維穩(wěn)態(tài)模型中的沉降系數(shù)。最后，通過統(tǒng)計方法，探求沉降系數(shù)和吸附系數(shù)二者之間的等效關(guān)系。

2.1 管道懸浮物運(yùn)移等參替代影響因素研究

2.1.1 管道懸浮物運(yùn)移室內(nèi)實驗

實驗裝置主要有帶有恒定水頭裝置的水箱、膠皮管道、閥門、水泵、針筒、穩(wěn)流箱、秒表、刻度尺等。實驗裝置圖如圖1所示。水箱帶有定水頭裝置，和閥門共同控制泥漿流出速度。管道長20 m，每隔2 m設(shè)置觀測點，觀測點處插入針筒方便采樣。

室內(nèi)實驗流程大致如下：預(yù)先配好一定濃度的泥漿懸浮液倒入水箱中，打開閥門讓泥漿流出，讓其在管道中運(yùn)移，同時不停地往水箱里加泥漿懸浮液使水箱內(nèi)水頭保持不變，運(yùn)移一定時間后，關(guān)掉出水口閥門和尾部閥門，迅速取樣。根據(jù)泥的質(zhì)量和泥漿懸浮液體積求得泥漿懸浮液濃度。實驗參數(shù)表如表1所示。

圖1 室內(nèi)實驗裝置示意圖

實驗編號123粒徑a/mm0.150.150.15管徑D/mm193876初始濃度C0/mg/L506050605766速度v/m/h430402402

2.1.2 管道懸浮物運(yùn)移數(shù)值模擬

根據(jù)實驗裝置及實驗參數(shù)設(shè)定，用FEFLOW軟件建立管道實驗?zāi)Ｐ?。模型示意圖如圖2所示，選擇管徑19 mm、管長20 m的模型，采用三角形網(wǎng)格剖分，共10 000個網(wǎng)格，設(shè)置與實驗對應(yīng)的10個觀測點(每隔2 m一個)。

圖2 管道數(shù)值模型示意圖

在管道左端設(shè)置與實驗相同定濃度，流速通過水頭和滲透系數(shù)控制，管道管徑和長度、模擬時間均與實驗一致。調(diào)試FEFLOW溶質(zhì)運(yùn)移模型吸附系數(shù)[10]，將FEFLOW模型計算結(jié)果與每個觀測點的懸浮物濃度觀測數(shù)據(jù)進(jìn)行擬合，得出10個與實驗觀測點濃度對應(yīng)的吸附系數(shù)。類似地，調(diào)試沉降系數(shù)，擬合一維河流穩(wěn)態(tài)模型計算結(jié)果與實驗觀測結(jié)果，得到對應(yīng)的沉降系數(shù)。

2.1.3 實驗結(jié)果

在盡量保持初始濃度、速度、粒徑相同條件下，實驗初步探究了管徑對沉降系數(shù)和吸附系數(shù)的影響，得出：管徑與沉降系數(shù)和吸附系數(shù)成反比，沉降系數(shù)和吸附系數(shù)與運(yùn)移距離無關(guān)。

圖3 實例條件下沉降系數(shù)與管徑的對應(yīng)關(guān)系圖

2.2 野外實例中管道懸浮物運(yùn)移等參替代關(guān)系構(gòu)建

實驗得出同種條件下吸附系數(shù)、沉降系數(shù)與距離無關(guān)，且得到三組不同管徑對應(yīng)的吸附系數(shù)和沉降系數(shù)。另外，根據(jù)巖溶地區(qū)另一個排泥庫發(fā)生的事故泥漿泄露資料求得的沉降系數(shù)為8.14×10-3(1/d)，根據(jù)FEFLOW模型模擬結(jié)果得到對應(yīng)吸附系數(shù)為0.08。而該排泥庫與論文研究的排泥庫相比主要是暗河管道大小不同，因此將三組室內(nèi)得到的沉降系數(shù)和吸附系數(shù)取結(jié)合另一排泥庫泄露時求得的對應(yīng)數(shù)據(jù)，統(tǒng)計得到沉降系數(shù)與管徑的對應(yīng)關(guān)系圖(見圖3)，以及吸附系數(shù)與管徑的對應(yīng)關(guān)系(見圖4)。

圖4 實例條件下吸附系數(shù)與管徑對應(yīng)關(guān)系

根據(jù)實際勘查資料，實例中地下巖溶管道約為5 m寬，外延圖3和圖4所得的關(guān)系，或代入圖中擬合公式，得實例沉降系數(shù)值為8.29×10-3(1/d)，對應(yīng)的吸附系數(shù)值為0.060 8。

3 巖溶管道等參替代數(shù)值模擬方法應(yīng)用實例

3.1 排泥庫工程概況

排泥庫位于廣西靖西縣，本區(qū)斷裂構(gòu)造發(fā)育，北西向斷層遍及全區(qū)，其中右江大斷裂F1規(guī)模最大，為控相區(qū)域性斷裂，把靖西～田東隆起帶錯移20多 km。根據(jù)野外巖溶水點調(diào)查及區(qū)域水文地質(zhì)資料分析，該區(qū)域地下暗河管道發(fā)育非常強(qiáng)烈，有龐凌(凌晚)地下河和鵝泉地下河兩個巖溶水補(bǔ)逕排系統(tǒng)，龐凌地下河與鵝泉地下河兩個地下河系統(tǒng)之間存在北西—南東向，次一級地下水分水嶺。

鵝泉地下河由1條主干支流和北面3條支流、南面1條支流構(gòu)成，補(bǔ)給面積為258 km2。地下河系沿巖溶管道、裂隙徑流，主干流程約21 km。 擬建排泥庫位于鵝泉地下河系統(tǒng)范圍內(nèi)，距鵝泉河干流約2.3 km。

圖5 模型3D示意圖

3.2 基于FEFLOW的懸浮物運(yùn)移數(shù)值模擬

正確可信的水流模型是溶質(zhì)運(yùn)移模型的前提和基礎(chǔ)，對水流模型進(jìn)行識別驗證。

3.2.1 水流模擬模型求解及識別驗證

根據(jù)研究區(qū)地質(zhì)及水文地質(zhì)勘察資料、示蹤試驗、鉆孔等資料分析污染物泥漿運(yùn)移的主要途徑，確定模擬計算區(qū)范圍為一相對獨立的水文地質(zhì)單元，西至F1斷裂為邊界，北至次一級地下水分水嶺，南邊界和東邊界為鵝泉地下暗河，面積約76 m2。模型示意圖如圖5所示。模型采用三角形單元網(wǎng)格剖分，研究區(qū)的總面積為76 km2，在地下暗河及支流處加密網(wǎng)格，

最小的網(wǎng)格面積為35 m2，最大為30 600 m2。在垂向上剖分為2層。研究區(qū)總計20 412個網(wǎng)格，16 218個節(jié)點。

根據(jù)巖性出露情況、巖層及地下暗河分布情況，將滲透系數(shù)(Kxx、Kyy)劃分為5個分區(qū)，巖性分區(qū)如圖6所示，對應(yīng)參數(shù)見表2。

圖6 滲透系數(shù)(Kxx、Kyy)分區(qū)圖

參數(shù)取值貯水系數(shù)/m/d0.001流入量/m/d0.0002孔隙度0.3滲透系數(shù)Kxx、Kyy/m/dI區(qū)0.50II區(qū)3.62III區(qū)2.45IV區(qū)7.54V區(qū)1000滲透系數(shù)Kzz/m/d0.01

模型識別驗證期從2014年3月至2015年3月，劃分為12個模擬時段，為與豐平枯三期的實測資料相對應(yīng)，時間步長為20～40d不等。應(yīng)用試估-校正法和單純形法識別得到模型參數(shù)(表)，所有時段水頭擬合誤差絕對值的均值在4.36 m以內(nèi)，考慮到巖溶地水文地質(zhì)條件的復(fù)雜性以及局部峰叢和溶洞對水流的影響，局部水位在現(xiàn)有條件下會有較大出入，但總體水位與實測情況吻合。

3.2.2 基于FEFLOW懸浮物運(yùn)移數(shù)值預(yù)測

污染物遷移模擬模型用于預(yù)測評價污染物在巖溶管道(地下暗河)的運(yùn)移情況，巖溶管道彌散度參數(shù)參考前人在巖溶管道進(jìn)行的相關(guān)野外試驗研究成果，縱向彌散度取值為400 m。污染泄露方式、濃度、等參替代模擬模型計算懸浮物運(yùn)移距離結(jié)果見表3。

3.3 懸浮物遷移預(yù)測結(jié)果分析

根據(jù)預(yù)測結(jié)果，基于FEFLOW的數(shù)值模擬結(jié)果總體略大于一維穩(wěn)態(tài)水質(zhì)模型解析法的結(jié)果，枯水期相差約500～800 m，豐水期相差約1 200～1 800 m左右。通過進(jìn)一步分析，解析法公式中，污染物濃度的影響因素只考慮初始濃度、流速、沉降系數(shù)三因素，不考慮水流的源匯項，且流速敏感性高于初始濃度的敏感性；而FEFLOW有限元法在考慮上述三項因素基礎(chǔ)上，考慮了河流與地下水之間的水力交換而產(chǎn)生的源匯項，因此懸浮物運(yùn)移距離稍大，與實際條件更為接近，亦更合理。

表3 泥漿懸浮物運(yùn)移距離計算結(jié)果一覽表 m

注：影響距離為泥漿濃度降至30 mg/L(根據(jù)鋁工業(yè)污染物排放標(biāo)準(zhǔn)GB 25465-2010水污染懸浮物限制為30 mg/L)的范圍。

4 結(jié)語

(1)FEFLOW模擬巖溶管道中懸浮類污染物遷移時，可以采用機(jī)理相似的吸附系數(shù)等效替代懸浮污染物運(yùn)移過程中的沉降系數(shù)。

(2)相比于常用的一維穩(wěn)態(tài)水質(zhì)模型，等參替代模型考慮了巖溶管道與周圍裂隙地下水之間的水力交換，模擬巖溶管道懸浮污染物運(yùn)移更為合理。

(3)影響溶質(zhì)運(yùn)移沉降作用的因素有很多，本文所探究的數(shù)值法模型采用的吸附系數(shù)和解析法參數(shù)沉降系數(shù)間的對應(yīng)關(guān)系僅是大量算例條件下得出，需要通過試驗和野外實際算例的進(jìn)一步研究和驗證。

[1]Shuster, E.T. and W.B. White, Seasonal fluctuations in the chemistry of lime-stone springs: A possible means for characterizing carbonate aquifers. Journal of Hydrology, 1971. 14(2): p. 93-128.

[2]Abusaada M, Sauter M. Studying the flow dynamics of a karst aquifer system with an equivalent porous medium model[J]. Groundwater, 2013.51(4): 641-650.

Fleury P, Plagnes V, Bakalowicz M. Modelling of the functioning of karst aquifers with a reservoir model: Application to Fontaine de Vaucluse (South of France)[J]. Journal of hydrology, 2007, 345(1): 38-49.

[3]Atkinson T C. Present and future directions in karst hydrogeology[J]. Annales de la Société géologique de Belgique, 1985.

[4]Quinlan J F, Ewers R O. Ground water flow in limestone terranes: Strategy rationale and procedure for reliable, efficient monitoring of ground water quality in karst areas[C]//Proceedings of the Fifth National Symposium and Exposition on Aquifer Restoration and Ground Water Monitoring May 21-24, 1985.The Fawcett Center, Columbus, Ohio. 1985. p 197-234, 12 fig, 100 ref. 1985.

[5]Scanlon, B.R., et al., Can we simulate regional groundwater flow in a karst system using equivalent porous media models Case study, Barton Springs Edwards aquifer, USA. Journal of Hydrology, 2003. 276(1-4): p. 137-158.

[6]薛禹群，吳吉春.地下水?dāng)?shù)值模擬在我國——回顧與展望. 水文地質(zhì)工程地質(zhì).1997. 24(4): 21-24.

[7]郝芳華, 李春暉, 趙彥偉. 流域水質(zhì)模型與模擬[J]. 北京師范大學(xué)出版社.2008.

[8]薛禹群，與春紅.地下水?dāng)?shù)值模擬. 科學(xué)出版社.2007.

[9]Karickhoff, S.W., D.S. Brown and T.A. Scott, Sorption of hydrophobic pollutants on natural sediments. Water research, 1979. 13(3): p. 241-248.

[10]張志強(qiáng)，等.基于示蹤試驗的巖溶管道及水力參數(shù)定量解析.人民長江.2015. 46(11): 80-83.

StudyonnumericalmodelofsuspendedsedimenttransportandequalparametersubstitutioninKarstpipelines——taking the prediction of suspended sediment transport in a row of mud banks in Guangxi as an example

ZHOUYu-chao1，ZHUXiao-bin1，WUJi-chun1，WUJian-feng1，ZENGXian-kui1，GUOFei2

(1.Key Laboratory of supergene geochemistry, Ministry of Education/ School of Earth Science and engineering, Nanjing University, Department of Water Science, Nanjing,Jiangsu,210023；2.Key Laboratory of virtual geographic environment, Ministry of education, Nanjing Normal University,Nanjing,Jiangsu,210046)

The migration of karst underground river in the area of suspension simulation is always difficult. In practical applications, due to the role of traditional settlement numerical model without suspension parameters of the reaction, usually only consider the generalization of suspended particles in the underground river in the river by using one-dimensional water quality prediction model, the lack of other methods for effective verification. This paper combines experiment and numerical simulation, laboratory experiments on the transport of suspended solids in pipes under different diameters and initial concentrations are designed, the feasibility of the equivalent substitute of the one-dimensional steady-state water quality model settlement coefficient in a numerical model with similar preliminary investigation mechanism is presented,coupled with mud in guangxi karst area library engineering build specific adsorption coefficient and subsidence coefficient under the condition of isoparametric alternative relations, finally the isoparametric alternative numerical model used in the karst area rivers suspended solids transport simulation, the numerical results show that the numerical model with equivalent parameters is more reasonable than one-dimensional steady-state water quality model. This study enriches the suspension of karst underground river transport simulation method.

Karst；one-dimensional steady water quality model；numerical simulation model；adsorption coefficient；sedimentation coefficient

P641.2

A

1004-1184(2017)05-0003-04

2017-06-17

國家自然科學(xué)基金項目(U1503282，41372235，41571386)；國家重點研發(fā)計劃課題(2016YFC0402807)

周宇超(1991-)，女，江蘇江陰人，在讀碩士研究生，主攻方向：地下水?dāng)?shù)值模擬及環(huán)境研究。

祝曉彬(1980-)，男，安徽東至人，副教授，主要從事水循環(huán)和地下水環(huán)境方面教學(xué)及研究工作。