楊 楊，趙良杰,2，蘇春田，夏日元

(1.中國地質(zhì)科學(xué)院巖溶地質(zhì)研究所/國土資源部廣西巖溶動(dòng)力學(xué)重點(diǎn)實(shí)驗(yàn)室,廣西 桂林 541004；2.中國地質(zhì)大學(xué)(北京)水資源與環(huán)境學(xué)院,北京 100083)

在我國西南地區(qū)，多重巖溶含水介質(zhì)常由基巖裂隙和巖溶管道構(gòu)成，水動(dòng)力場(chǎng)非常復(fù)雜[1]，傳統(tǒng)的地下水?dāng)?shù)值模擬難以刻畫巖溶地下水流特征[2-3]。Shoemaker等在MODFLOW-2005開源程序的基礎(chǔ)上增加管道流程序(Conduit Flow Process，CFP)，分別應(yīng)用Hagen-Poiseuille方程和Darcy-Weisbach方程描述巖溶管道介質(zhì)中層流和紊流運(yùn)動(dòng)[4-5]。CFP能夠刻畫巖溶管道特征和快速的優(yōu)先流，適用于復(fù)雜的多重巖溶含水介質(zhì)。許多學(xué)者運(yùn)用理想模型對(duì)CFP性能進(jìn)行測(cè)試：Hill等[6]通過對(duì)比傳統(tǒng)等效多孔介質(zhì)模型和管道流CFP模型(雙重介質(zhì))，基于觀測(cè)和模擬流量結(jié)果比較，得出的管道流程序模擬更加精確；Gallegos等[7]基于物理沙箱模型和Woodville巖溶區(qū)觀測(cè)數(shù)據(jù)，運(yùn)用CFP模擬泉水流量，得出管道流程序能夠模擬巖溶泉流量動(dòng)態(tài)特征；Giese等[8]應(yīng)用CFP模型研究管道中水流的紊流運(yùn)動(dòng)特征；常勇[9]基于MODFLOW-CFP模型模擬巖溶含水層內(nèi)管道-裂隙二元結(jié)構(gòu)水文過程，并提出一種新的水箱——CFP組合模型；焦友軍等10]采用CFP模型模擬巖溶含水系統(tǒng)在暴雨期的響應(yīng)過程。

CFP模型也可用于研究溶質(zhì)運(yùn)移特征。Faulkner等[11]利用實(shí)驗(yàn)室物理模型和數(shù)值模型研究基巖裂隙和巖溶管道雙重介質(zhì)溶質(zhì)運(yùn)移特征和交換規(guī)律，將巖溶含水層分為兩個(gè)區(qū)域(基巖裂隙和巖溶管道區(qū))，地下水在基巖裂隙區(qū)符合達(dá)西定律，而在巖溶管道內(nèi)滿足Navier-Stokes方程，數(shù)值模擬的水頭分布、流量和溶質(zhì)運(yùn)移特征與物理實(shí)驗(yàn)結(jié)果一致。Ghasemizadeh等[12]對(duì)巖溶區(qū)等效介質(zhì)、裂隙介質(zhì)及管道介質(zhì)水流模型和溶質(zhì)運(yùn)移模型進(jìn)行總結(jié)。等效介質(zhì)和裂隙介質(zhì)模型所需數(shù)據(jù)少，概化巖溶介質(zhì)內(nèi)部結(jié)構(gòu)特征建模方便，但不能刻畫水量和水質(zhì)的空間分布特征；等效介質(zhì)模型不適應(yīng)于巖溶區(qū)，而在等效介質(zhì)基礎(chǔ)上耦合分布式管道(CFP)模型能夠模擬管道-裂隙介質(zhì)水流系統(tǒng)，但對(duì)模型數(shù)據(jù)和參數(shù)需求較高。模塊化三維溶質(zhì)遷移模型(The modular three-dimentional transport model,MT3DMS)，是國際上應(yīng)用最廣的地下水系統(tǒng)污染物遷移模擬程序[13-15]?；趲r溶管道流程序研究巖溶管道裂隙水流運(yùn)動(dòng)較為成熟[16-18]，但結(jié)合污染物在管道中的運(yùn)移特征研究相對(duì)較少。

本文的目的是以概念模型為例，通過耦合CFP和MT3DMS研究巖溶管道流溶質(zhì)運(yùn)移特征，并討論管道參數(shù)如管道直徑、水力梯度及水流交換系數(shù)對(duì)溶質(zhì)運(yùn)移模擬結(jié)果的影響，以期更好地理解實(shí)際巖溶區(qū)溶質(zhì)運(yùn)移現(xiàn)象。

1 研究方法

1.1 管道流模型(CFP)

CFP模型是在孔隙介質(zhì)層流模型(MODFLOW-2005)基礎(chǔ)上耦合離散的管道網(wǎng)絡(luò),符合西南巖溶區(qū)雙重介質(zhì)水流系統(tǒng)[19-20](圖1)。

圖1 CFP原理示意圖Fig.1 Schematic diagram showing the principle of CFP

孔隙介質(zhì)層流模型基于單元中心有限差分法,層流模型地下水在三維空間的流動(dòng)符合達(dá)西定律，公式如下：

(1)

式中：Kxx、Kyy、Kzz——滲透系數(shù)在x、y、z方向上的分量/(m·d-1)；

h——水頭/m；

W——單位體積流量，代表流進(jìn)匯和來自源的水量/d-1；

Ss——孔隙介質(zhì)的貯水率/m-1；

t——時(shí)間/d。

式(1)假設(shè)地下水流無黏性，密度不變，不可壓縮且屬于達(dá)西流。滲透系數(shù)的主軸方向與坐標(biāo)軸方向一致。加上相應(yīng)的初始條件和邊界條件,構(gòu)成地下水流孔隙介質(zhì)數(shù)學(xué)模型。MODFLOW-2005提供多種數(shù)值方法求解線性及非線性給定水頭邊界、給定流量邊界及第三類邊界條件。

管道系統(tǒng)是由節(jié)點(diǎn)和圓柱形管道組成,可參考Shoemaker等發(fā)表的CFP使用手冊(cè)。根據(jù)體積守恒定律,每一個(gè)節(jié)點(diǎn)水流狀態(tài)可應(yīng)用Kirchhoff定律描述(式2)。當(dāng)管道水流處于層流狀態(tài)時(shí),應(yīng)用Hagen-Poiseuille(式3)描述；管道水流處于紊流狀態(tài)時(shí),應(yīng)用Darcy-Weisbach(式4)描述，CFP應(yīng)用式(5)計(jì)算管道與巖溶基巖水流交換。

(2)

(3)

(4)

Qex=αex(hc-hm)

(5)

式中：Qip——從管道i到管道n的流量/(m3·d-1)；

Qss——所有節(jié)點(diǎn)源匯項(xiàng)流量總和/(m3·d-1)；

d——管道直徑/m；

g——重力加速度/(m·d-2)；

Δhc——管道水頭損失/m；

υ——水動(dòng)力黏滯系數(shù)/(m2·d-1)；

kc——管道平均粗糙度，與管道壁形態(tài)有關(guān)；

Qex——管道與基巖交換水量/(m3·d-1)，值為負(fù)表示水流從基巖流向巖溶管道，值為正表示水流從管道流向基巖含水層；

αex——管道滲透系數(shù)/(m2·d-1)；

hc——管道內(nèi)水頭/m；

hm——基巖水頭/m。

1.2 模塊化三維溶質(zhì)遷移模型(MT3DMS)

MT3DMS是國際上應(yīng)用最廣的地下水系統(tǒng)污染物遷移模擬程序(對(duì)流、彌散及生化反應(yīng)等),該模型能夠靈活有效地處理各種源匯項(xiàng)和邊界條件，適用于承壓、無壓及承壓-無壓含水層中的污染物運(yùn)移問題。偏微分式(6)用于描述污染物k在三維非穩(wěn)定流系統(tǒng)中的溶質(zhì)運(yùn)移特性。

于學(xué)生現(xiàn)有認(rèn)知程度上對(duì)學(xué)生科學(xué)素養(yǎng)進(jìn)行進(jìn)一步提升為高中物理教學(xué)重要教學(xué)目標(biāo).因高中生即將面臨高考，而物理學(xué)習(xí)難度同其他學(xué)科相比學(xué)習(xí)難度較大，讓學(xué)生對(duì)物理教學(xué)予以正確理解可有效幫助學(xué)生構(gòu)建系統(tǒng)性較強(qiáng)的物理思維體系.高中學(xué)生物理核心素養(yǎng)是將學(xué)生物理知識(shí)、物理學(xué)習(xí)態(tài)度、物理學(xué)習(xí)方式及學(xué)生物理操作能力為基礎(chǔ)所提出的學(xué)生核心素養(yǎng)之一.核心素養(yǎng)不僅包括物理理論知識(shí)，同時(shí)也包括學(xué)生于所開展的物理實(shí)踐過程，為高中學(xué)生提高物理綜合能力的科學(xué)導(dǎo)向.由上述內(nèi)容可知，培養(yǎng)高中學(xué)生科學(xué)價(jià)值觀及強(qiáng)化高中學(xué)生社會(huì)認(rèn)知為高中物理核心素養(yǎng)重要價(jià)值體現(xiàn).

(6)

式中：R——阻滯因子，無量綱；

θ——含水介質(zhì)孔隙度，無量綱；

Ck——水中的污染物k的濃度/(kg·m-3)；

t——時(shí)間/d；

xi、xj——空間坐標(biāo)/m；

Dij——水動(dòng)力彌散系數(shù)張量/(m2·d-1)；

vi——地下水滲透流速/(m·d-1)；

qs——源(正值)或匯(負(fù)值)的單位流量/d-1；

∑Rn——化學(xué)反應(yīng)項(xiàng)/(kg·m-3·d-1)。

MT3DMS采用模塊化的結(jié)構(gòu)子程序包來實(shí)現(xiàn)不同功能。不包括求解地下水流方程的子程序，通過任何塊中心有限差分水流模型的輸出結(jié)果來獲得地下水位資料，如MODFLOW等。其中子程序包主要包括基本運(yùn)移子程序包(BTN)、水流模型接口包(FMI)、對(duì)流子程序包(ADV)、彌散子程序包(DSP)、源匯項(xiàng)子程序包(SSM)、化學(xué)反應(yīng)子程序包(RCT)、GCG迭代求解程序包(GCG)及實(shí)用工具子程序包(UTL)。

2 概念模型

管道流CFP中概念模型分為兩層，含水層由基巖網(wǎng)格和管道組成。第一層表示潛水含水層，第二層表示承壓含水層，管道設(shè)置在第二層，每層厚20 m(圖2)。網(wǎng)格大小10 m×10 m，共計(jì)47列，31行。水平滲透系數(shù)為1 m/h,垂向各向異性比為10,潛水含水層給水度為0.1,承壓含水層儲(chǔ)水系數(shù)為0.000 1。模型第1列和第47列為給定水頭邊界，其余為零通量邊界。其中管道P1至管道P4，管道直徑分別為0.2 m、0.3 m、0.4 m和0.5 m，彎曲度為1，管道粗糙度為0.001，下臨界雷諾數(shù)為2 000，上臨界雷諾數(shù)為4 000。節(jié)點(diǎn)N1至節(jié)點(diǎn)N5滲透系數(shù)αex為0.1 m2/h。

圖2 管道流CFP概念模型Fig.2 Conceptual model for CFP

應(yīng)力期分為8個(gè)，暴雨期前10 h平水期為第一應(yīng)力期，第二至第七應(yīng)力期暴雨期間隔4 h，暴雨后持續(xù)模擬686 h，共計(jì)30 d，落水洞處N1節(jié)點(diǎn)為集中降雨補(bǔ)給，其余為地表降水面狀補(bǔ)給。在沒有降雨情況下，第二層含水介質(zhì)等水位線圖如圖3所示，其中第1、47列為給定水頭邊界，水流自左向右流動(dòng)，等水位線在管道處有突變。由于水力梯度大，管道處于承壓狀態(tài)，雷諾數(shù)最小為37 047，水流為紊流狀態(tài)。

圖3 巖溶含水層第二層等水位線圖Fig.3 Contours of groundwater level in the 2th layer of the karst aquifer

從表1可以看出，在模型末期，節(jié)點(diǎn)N1、N2、N3、N4處管道水頭低于周圍基巖含水層水頭(N5為給定水頭)，水流從周圍含水層向管道徑流，管道為承壓狀態(tài)，其中節(jié)點(diǎn)N1接受含水層補(bǔ)給為187.10 m3/h，并流向管道P1；節(jié)點(diǎn)N2接受含水層補(bǔ)給296.29 m3/h，管道P1補(bǔ)給187.10 m3/h，并流向管道P2；節(jié)點(diǎn)N3接受含水層補(bǔ)給224.39 m3/h，管道P2補(bǔ)給483.39 m3/h，并流向管道P3；節(jié)點(diǎn)N4接受含水層補(bǔ)給124.24 m3/h，管道P3補(bǔ)給707.78 m3/h，并流向管道P4，節(jié)點(diǎn)N5流出量為832.02 m3/h，滿足守恒定律(式5)。在此水流模型基礎(chǔ)上運(yùn)行MT3DMS模型，為模擬巖溶區(qū)示蹤試驗(yàn)，僅在落水洞處第一應(yīng)力期示蹤劑濃度為100 mg/L，其余初始濃度為0。不考慮化學(xué)反應(yīng)情況下，選擇對(duì)流、彌散及源匯項(xiàng)子程序包?？紫抖葹?.1，縱向彌散度為100 m，橫向彌散度為10 m，垂向彌散度為1 m。Field計(jì)算出巖溶地區(qū)彌散系數(shù)介于0.08～1 m2/s[21]，而在實(shí)際工作中，通過示蹤試驗(yàn)計(jì)算巖溶管道內(nèi)流速介于0.002～1.121 m/s(計(jì)算管道平均直徑和出口流量)，彌散度介于1～500 m之間，因此彌散度設(shè)置為100 m是合理的。

表1 模擬末期水流交換特征Table 1 Characteristics of water exchange between conduit and matrix

分別對(duì)比模擬N2至出口N5處濃度運(yùn)移曲線(圖4)，N2、N3、N4、N5處距離投放點(diǎn)距離分別為100,200,300,400 m。從圖4中可以看出，管道中溶質(zhì)運(yùn)移存在明顯拖尾現(xiàn)象，且距離越長拖尾越明顯；距離越短峰值濃度越大，到達(dá)峰值的時(shí)間越短，符合實(shí)際巖溶區(qū)示蹤試驗(yàn)穿透曲線。

圖4 不同距離接收點(diǎn)濃度曲線Fig.4 Concentration curves at different distances

3 討論

為研究不同水文地質(zhì)參數(shù)對(duì)溶質(zhì)運(yùn)移的影響，選取管道P3中點(diǎn)距離投放點(diǎn)250 m處為研究對(duì)象，探討孔隙度、彌散系數(shù)、水力梯度、管道直徑和管道滲透系數(shù)對(duì)溶質(zhì)穿透曲線的影響。由于巖溶區(qū)濃度穿透曲線存在拖尾現(xiàn)象，為評(píng)估曲線的拖尾情況，引入偏態(tài)系數(shù)SK，該參數(shù)是研究曲線不對(duì)稱程度的統(tǒng)計(jì)參數(shù)(式7)。當(dāng)偏態(tài)系數(shù)為0時(shí)，表示濃度穿透曲線完全對(duì)稱分布；當(dāng)偏態(tài)系數(shù)大于0時(shí)，表示濃度穿透曲線為正偏態(tài)；當(dāng)偏態(tài)系數(shù)小于0時(shí)，表示濃度穿透曲線為負(fù)偏態(tài)。

(7)

式中：SK——偏態(tài)系數(shù)；

Ci——第i個(gè)濃度值/(kg·m-3)；

σ——均方差；

N——個(gè)數(shù)。

3.1 含水層孔隙度對(duì)溶質(zhì)運(yùn)移的影響

保持其它參數(shù)不變，選取孔隙度分別為0.250、0.100，0.075時(shí)，對(duì)比不同孔隙度接收點(diǎn)濃度曲線變化(圖5)。從圖5可以看出，當(dāng)θ為0.25和0.1時(shí)，濃度穿透曲線輕微不對(duì)稱，而θ當(dāng)為0.075時(shí)，濃度曲線顯著不對(duì)稱，且拖尾較長。隨著孔隙度的增大，濃度曲線峰值依次為12.51，30.68，40.53 mg/L，峰值到達(dá)時(shí)間依次為86，38，30 h，曲線歷時(shí)時(shí)間依次為710，484，388 h，偏態(tài)系數(shù)分別為1.58，1.25，0.98?？紫抖仍叫?，濃度曲線峰值越大，峰值到達(dá)時(shí)間越快，曲線歷時(shí)時(shí)間越短，濃度穿透曲線越對(duì)稱。

圖5 不同孔隙度接收點(diǎn)濃度曲線Fig.5 Concentration curves with different porosities

3.2 不同水力梯度對(duì)溶質(zhì)運(yùn)移的影響

保持其它參數(shù)不變，選取水力梯度i分別為0.01，0.04，0.10時(shí)，對(duì)比不同水力梯度接收點(diǎn)濃度曲線變化(圖6)。隨著水力梯度增大，濃度曲線峰值依次為4.55，30.68，82.52 mg/L，峰值到達(dá)時(shí)間依次為134，38，18 h，曲線歷時(shí)時(shí)間依次為720，484，220 h，偏態(tài)系數(shù)分別為0.74，1.25，0.52。水力梯度越大，濃度曲線峰值越大，峰值到達(dá)時(shí)間越快，曲線歷時(shí)時(shí)間越短。

圖6 不同水力梯度接收點(diǎn)濃度曲線Fig.6 Curves with different hydraulic gradients

3.3 管道直徑對(duì)溶質(zhì)運(yùn)移的影響

保持其它參數(shù)不變，選取管道直徑d分別為0.01，0.10，0.50時(shí)，對(duì)比不同管道直徑接收點(diǎn)濃度曲線變化(圖7)。隨著管道直徑增大，濃度曲線峰值依次為16.09，27.56，29.32 mg/L，峰值到達(dá)時(shí)間依次為38，34，38 h，曲線歷時(shí)時(shí)間依次為412，460，508 h，偏態(tài)系數(shù)分別為1.02，1.14，1.30。管道直徑越大，濃度曲線峰值越大，峰值到達(dá)時(shí)間基本相同。

圖7 不同管道直徑接收點(diǎn)濃度曲線Fig.7 Curves with different conduit diameters

3.4 管道滲透系數(shù)αex對(duì)溶質(zhì)運(yùn)移的影響

保持其它參數(shù)不變，選取管道滲透系數(shù)αex分別為0.01，0.10，1.00時(shí)，對(duì)比不同管道滲透系數(shù)接收點(diǎn)濃度曲線變化(圖8)。隨著管道滲透系數(shù)增大，濃度曲線峰值依次為16.59，30.68，30.02 mg/L，峰值到達(dá)時(shí)間依次為38，38，38 h，曲線歷時(shí)時(shí)間依次為412，487，604 h，偏態(tài)系數(shù)分別為1.04，1.25，1.37。管道滲透系數(shù)越大，濃度曲線峰值越大，曲線越不對(duì)稱，峰值到達(dá)時(shí)間相同。

圖8 不同管道滲透系數(shù)接收點(diǎn)濃度曲線Fig.8 Curves with different pipe conductances

4 結(jié)論

(1)通過概念模型算例，研究不同水文地質(zhì)參數(shù)對(duì)溶質(zhì)運(yùn)移的影響，發(fā)現(xiàn)隨著水力梯度、管道直徑及管道滲透系數(shù)的增大，孔隙度減小，濃度曲線峰值越大，峰值到達(dá)時(shí)間越快，濃度穿透曲線越對(duì)稱，為將來研究實(shí)際巖溶區(qū)溶質(zhì)運(yùn)移現(xiàn)象，如示蹤試驗(yàn)穿透曲線、尾礦庫污染物在地下河出口的變化規(guī)律等提供科學(xué)基礎(chǔ)。

(2)管道流CFP模型能夠刻畫巖溶管道與基巖裂隙水流交換特征，耦合MT3DMS溶質(zhì)運(yùn)移模型能夠模擬穿透曲線的拖尾現(xiàn)象，符合實(shí)際巖溶區(qū)特征，為分析巖溶管道流溶質(zhì)運(yùn)移特征提供了一種有效手段。

(3)實(shí)際巖溶區(qū)溶質(zhì)運(yùn)移濃度穿透曲線經(jīng)常存在多峰曲線，管道中存在溶潭、水流陡坎等不連續(xù)現(xiàn)象，而CFP模型中僅用平直圓管概化巖溶管道特征，與實(shí)際巖溶管道形態(tài)有一定出入；精確刻畫巖溶管道形態(tài)特征是進(jìn)一步研究的目標(biāo)。