郝奇琛，崔偉哲，黃林顯

(1. 中國(guó)地質(zhì)科學(xué)院水文地質(zhì)環(huán)境地質(zhì)研究所，河北石家莊050061； 2. 河北地質(zhì)大學(xué)水資源與環(huán)境學(xué)院, 河北石家莊050031；3. 濟(jì)南大學(xué)水利與環(huán)境學(xué)院，山東濟(jì)南250022)

盆地地下水流系統(tǒng)的形成與演變受多種因素的控制和影響[1-2]。根據(jù)地下水流系統(tǒng)理論，地形勢(shì)能是地下水流系統(tǒng)形成的主要驅(qū)動(dòng)因子，即地形控制論。事實(shí)上，除地形外，密度也是驅(qū)動(dòng)地下水流動(dòng)的因子之一，并且逐漸引起學(xué)者們的關(guān)注[3-5]。

海水入侵現(xiàn)象是密度影響水流運(yùn)動(dòng)最為典型的例子，密度較大的咸水會(huì)侵入到海岸帶含水層底部[6-7]。海水與淡水的密度差只有2.5%，干旱盆地中部由于蒸發(fā)濃縮作用，因此地下鹵水與淡水的密度差最大可達(dá)20%～30%，加之沉積盆地中部多分布著低滲透介質(zhì)，變密度影響下的地下水運(yùn)動(dòng)變得更為復(fù)雜[8-9]。Zech等[10]對(duì)德國(guó)某小型不規(guī)則沉積盆地地下水流系統(tǒng)進(jìn)行了模擬研究，發(fā)現(xiàn)高密度對(duì)區(qū)域水流系統(tǒng)影響較大，有明顯的削弱作用，但未詳細(xì)闡明產(chǎn)生這一現(xiàn)象的原因。

本文中主要采用平均等價(jià)淡水水頭差分析方法，構(gòu)建水流驅(qū)動(dòng)力削弱程度的定量表征模型，基于典型剖面的變密度流數(shù)值模擬結(jié)果，闡述密度變化對(duì)區(qū)域水流系統(tǒng)地下水徑流路徑的影響，以及對(duì)區(qū)域水流驅(qū)動(dòng)力的削弱機(jī)制。

1 研究區(qū)概況

研究區(qū)位于我國(guó)西部的柴達(dá)木盆地，是典型的內(nèi)陸封閉盆地。氣候?qū)賰?nèi)陸極端干旱高寒氣候，盆地周圍高山環(huán)繞，暖濕氣流難以進(jìn)入，所以降水稀少，蒸發(fā)量大。據(jù)當(dāng)?shù)貧庀笳径嗄暧^測(cè)資料分析，盆地多年平均降水量為16.09～189.73 mm，平原區(qū)多年平均蒸發(fā)量為1 973.62～3 183.04 mm。區(qū)內(nèi)河流主要分布有那陵格勒河、格爾木河、香日德河、巴音河、塔塔棱河、柴達(dá)木河、大哈勒騰河、小哈勒騰河等。

在盆地的形成過(guò)程中，大量的松散沉積物堆積在盆地內(nèi)，沉積物的孔隙結(jié)構(gòu)為地下水的儲(chǔ)存提供了良好的儲(chǔ)水空間。受水動(dòng)力條件影響，顆粒較粗的砂卵礫石堆積在山前的沖洪積扇地帶，形成一層或多層厚度約幾十米的含水介質(zhì)，有的厚度高達(dá)數(shù)百米。山前沖洪積扇含水層滲透性好，富水性強(qiáng)，地下水徑流速度較快。從沖洪積扇中后部到前緣的溢出帶，地層巖性由單一的砂卵礫石層逐漸變?yōu)槎鄬由芭c亞砂土互層，由單層結(jié)構(gòu)變?yōu)槎鄬咏Y(jié)構(gòu)，單個(gè)含水層厚度變薄，滲透性能減弱，地下水更新速率變慢，地下水類型由潛水過(guò)渡為上部潛水、下部承壓水。盆地中部沉積物厚度變大，但有效含水層厚度或單個(gè)含水層厚度變小，一般為幾米或者幾十米，富水性相對(duì)較差，含水介質(zhì)主要為粉細(xì)砂，含水層之間分布多層亞砂土及黏土等弱透水層，滲透性能變得更差。

從沖洪積扇到溢出帶，再到湖積平原，地下水補(bǔ)給、徑流、排泄條件發(fā)生規(guī)律性的變化，形成多級(jí)地下流系統(tǒng)。從補(bǔ)給區(qū)到排泄區(qū)地下水水化學(xué)分帶特征明顯，水化學(xué)組分及密度差異極大，是研究水化學(xué)組分變化與地下水流系統(tǒng)相互作用的理想?yún)^(qū)域[11]。

2 研究方法

2.1 數(shù)值模型設(shè)置

根據(jù)盆地水文地質(zhì)條件特征，將理想模型的剖面長(zhǎng)度設(shè)定為60 km，盆地中部最大深度為1 500 m，盆地邊緣最小厚度為200 m。為了便于分析，理想模型設(shè)置為對(duì)稱結(jié)構(gòu)(見圖1)。模型采用不規(guī)則六面體的網(wǎng)格剖分方式，網(wǎng)格的水平間距為1 000 m，共剖分50層，垂向間距由地表向下逐漸增大，最小約為0.10 m，最大約為100 m，剖分的網(wǎng)格數(shù)為3 000(見圖1)。

從盆地邊緣到盆地中部，大致可概化為3個(gè)地貌分帶，即沖洪積扇、溢出帶、湖積平原。巖性組合分別為中砂、中砂-細(xì)砂、細(xì)砂-亞砂土。模型中使用的水文地質(zhì)參數(shù)根據(jù)巖性分區(qū)賦值；溶質(zhì)運(yùn)移參數(shù)及水文地球化學(xué)參數(shù)，均采用經(jīng)驗(yàn)值并根據(jù)模型模擬結(jié)果進(jìn)行合理校正。表1所示為識(shí)別后的水文地質(zhì)參數(shù)，垂向滲透系數(shù)取水平滲透系數(shù)的1/10，擴(kuò)散系數(shù)取值為3.921×10-6m2/s。

模擬區(qū)地下水唯一的補(bǔ)給來(lái)源為河流入滲補(bǔ)給，補(bǔ)給量大小根據(jù)內(nèi)陸盆地河流入滲量大致確定。主要排泄項(xiàng)為泉、河流、潛水蒸發(fā)。理想模型主要針對(duì)一般規(guī)律的研究，忽略人為開采的影響。地下水系統(tǒng)水溶組分的補(bǔ)給項(xiàng)包括隨河流入滲至地下水的水溶組分、地下水溶濾含水介質(zhì)的組分。主要的排泄項(xiàng)包括隨泉水及河水流出的組分、發(fā)生礦物沉淀的組分。

(a)典型剖面平面位置圖

(b)典型剖面巖性分區(qū)示意圖

(c)典型剖面網(wǎng)格剖分示意圖圖中的地理數(shù)據(jù)由國(guó)家基礎(chǔ)地理信息中心網(wǎng)站下載(http://www.webmap.cn/mapDataAction.do?method=forw&resType=5&storeId=2&storeName=%E5%9B%BD%E5%AE%B6%E5%9F%BA%E7%A1%80%E5%9C%B0%E7%90%86%E4%BF%A1%E6%81%AF%E4%B8%AD%E5%BF%83&fileId=BA420C422A254198BAA5ABAB9CAAFBC1)；地勢(shì)數(shù)據(jù)從地理空間數(shù)據(jù)云網(wǎng)站下載(http://www.gscloud.cn/search)，利用ArcGIS10.3軟件數(shù)字化處理后得到。圖1 剖面地質(zhì)結(jié)構(gòu)概化及網(wǎng)格剖分示意圖

表1 識(shí)別后的水文地質(zhì)參數(shù)

表2 組分初始濃度值

水鹽的富集過(guò)程是一個(gè)長(zhǎng)期的地質(zhì)演化過(guò)程，因此將模型的應(yīng)力期定為100 ka。最少模擬步長(zhǎng)設(shè)置為1 s，最大為10 ka。本文中的模擬工作采用的Toughreact程序具有自動(dòng)調(diào)節(jié)時(shí)間步長(zhǎng)功能，會(huì)根據(jù)收斂難易程度來(lái)自動(dòng)調(diào)整時(shí)間步長(zhǎng)的大小[12]。

2.2 水流驅(qū)動(dòng)力削弱程度定量表征方法

驅(qū)動(dòng)力反映了地下水從補(bǔ)給區(qū)到排泄區(qū)勢(shì)能降低的趨勢(shì)，是一種自組織行為。地下水動(dòng)力學(xué)以水頭表示地下水總勢(shì)能狀態(tài)?；谶@一經(jīng)典方法，本文中同樣以水頭差來(lái)定量表征水流驅(qū)動(dòng)力的大小。

圖2(a)所示為淡水條件下山前補(bǔ)給區(qū)到盆地中心排泄區(qū)的區(qū)域水流概念模型，其中H1、H2分別為補(bǔ)給區(qū)和排泄區(qū)靜水條件下的水頭值，底部為基準(zhǔn)面，ΔH為兩者之間的差值，可以表征兩點(diǎn)之間地形驅(qū)動(dòng)力的大小。圖2(b)所示為排泄區(qū)為鹵水條件下的概念模型。由于盆地排泄區(qū)地下水密度增大，H1與H2的差值不再適合表示驅(qū)動(dòng)力的大小，因此需要將含鹽水頭值H2換算為等價(jià)淡水水頭值Hf，用H1與Hf的差值表示地形驅(qū)動(dòng)力的大小。

等價(jià)淡水水頭與含鹽水頭的換算公式為

(1)

式中:ρs、ρf分別為鹽水和淡水密度，kg/m3；Z為參考基準(zhǔn)面，m。

參照Hamann等[9]對(duì)平均密度差的表達(dá)方法，本文中擬采用平均等價(jià)淡水水頭差(average head contrast，AHC)δAHC表示密度增大導(dǎo)致的排泄區(qū)水頭增加值，即

(a)淡水條件

(b)鹵水條件H1—補(bǔ)給區(qū)靜水條件下的水頭值； H2—排泄區(qū)靜水條件下的水頭值； Hf—等價(jià)淡水水頭值； ΔH—H1與H2的差值； ΔHf—Hf與H2的差值； W—盆地排泄區(qū)寬度。圖2 盆地補(bǔ)給區(qū)到排泄區(qū)區(qū)域水流概念模型

(2)

式中：W、L分別為盆地排泄區(qū)寬度和深度，m； ΔHf為等價(jià)淡水水頭差，m，是鹽水的等價(jià)淡水水頭(圖2(b)中Hf)與含鹽水頭(圖2(b)中H2)的差值，隨時(shí)間和空間發(fā)生變化。

同理，設(shè)淡水條件下補(bǔ)給區(qū)的平均水頭值為H1a，排泄區(qū)的平均水頭值為H2a，兩者差值即可表征水流驅(qū)動(dòng)力的大小。排泄區(qū)密度造成的水頭增加值除以H1a與H2a的差值，即可表征水流驅(qū)動(dòng)力的削弱程度E，即

E=[δAHC/(H1a-H2a)]×100%。

(3)

借助以上定量表征方法，根據(jù)理想剖面變密度流模型的模擬結(jié)果，可以定量分析排泄區(qū)密度增大對(duì)水流驅(qū)動(dòng)力的削弱程度。

3 結(jié)果與討論

3.1 地下水流系統(tǒng)劃分

地下水流系統(tǒng)是研究水質(zhì)(鹽量、熱量)時(shí)空演變的理想框架與工具。根據(jù)地下水流系統(tǒng)理論，流域盆地可發(fā)育嵌套式的多級(jí)水流系統(tǒng)，包括局部、中間及區(qū)域水流系統(tǒng)，地下水流以不同級(jí)次方式有序運(yùn)移，水量、鹽分、熱量發(fā)生有規(guī)律的時(shí)空演變，呈現(xiàn)出時(shí)空有序的結(jié)構(gòu)。地下水流系統(tǒng)的劃分對(duì)研究地下水循環(huán)具有重要意義。

基于Tóth[13]提出的以2條并不總是相鄰的流線作為不同流動(dòng)系統(tǒng)的邊界的方法，利用Tough-react程序的模擬結(jié)果，生成水頭及地下水流線圖(見圖3)，再結(jié)合地形、沉積條件、水文地質(zhì)參數(shù)劃分出3個(gè)級(jí)別地下水流系統(tǒng)，即局部水流系統(tǒng)、中間水流系統(tǒng)、區(qū)域水流系統(tǒng)，如圖4所示。在山前接受補(bǔ)給、在溢出帶以泉的形式排泄的地下水可劃分為一個(gè)獨(dú)立的水流系統(tǒng)，稱為山前局部水流系統(tǒng)，記為Ⅰ1。溢出帶之后，部分泉集河水入滲到地下，并以蒸發(fā)的形式排泄，也可劃分為一個(gè)水流系統(tǒng)，稱為溢出帶局部水流系統(tǒng)，記為Ⅰ2。部分山前補(bǔ)給的地下水跨越2個(gè)局部水流系統(tǒng)后以潛水蒸發(fā)或構(gòu)造泉的形式排泄，形成中間水流系統(tǒng)，記為Ⅱ。從山前接受補(bǔ)給，跨越中間水流系統(tǒng)，到盆地中部湖積平原區(qū)以蒸發(fā)的形式排泄，最終形成區(qū)域水流系統(tǒng)，記為Ⅲ。

圖3 模擬的水頭及流線分布

圖4 盆地地下水流系統(tǒng)劃分

本文中重點(diǎn)討論地下水密度增加對(duì)區(qū)域水流系統(tǒng)的影響。經(jīng)過(guò)初步計(jì)算，發(fā)現(xiàn)地下水密度增加對(duì)水流驅(qū)動(dòng)力的削弱程度在區(qū)域水流系統(tǒng)的上層和下層差異較大，因此又將區(qū)域水流系統(tǒng)細(xì)分為上層區(qū)域水流系統(tǒng)和下層區(qū)域水流系統(tǒng)。

3.2 地下水密度的變化

3.2.1 剖面密度變化

圖5所示為模擬的不同時(shí)期地下水密度的變化。在模型的初始時(shí)刻，地下水處于淡水狀態(tài)。蒸發(fā)濃縮是干旱盆地主要的水文地球化學(xué)作用，受蒸發(fā)作用的影響，盆地中部排泄區(qū)鹽分富集，地下水密度逐漸增大。鹽分由盆地中心逐漸向下、向兩側(cè)擴(kuò)散。模型運(yùn)行至10 ka時(shí)，盆地中部地下水密度最大值為1.01×103kg/m3，溶解性總固體(TDS)質(zhì)量濃度最大值為9 g/L，局部變?yōu)橄趟臀⑾趟?，此時(shí)盆地中部地下水與淡水密度差都在1%以下。模型運(yùn)行至30 ka，盆地中部地下水密度最大值為1.05×103kg/m3，TDS質(zhì)量濃度最大值為47 g/L，局部已經(jīng)變?yōu)辂}水。模型運(yùn)行至50 ka，地下水密度最大值為1.08×103kg/m3，TDS質(zhì)量濃度最大值為65 g/L，此時(shí)盆地中部大部分地區(qū)已經(jīng)大于鹽水的密度，部分地區(qū)開始出現(xiàn)鹵水。模型運(yùn)行至70～100 ka，盆地中部出現(xiàn)大面積鹵水，地下水密度最大值為1.11×103kg/m3，TDS質(zhì)量濃度最大值為98 g/L，此時(shí)盆地中部地下水與淡水密度差最大為12%。

(a)初始時(shí)刻(b)模擬時(shí)間為10 ka(c)模擬時(shí)間為30 ka(d)模擬時(shí)間為50 ka(e)模擬時(shí)間為70 ka(f)模擬時(shí)間為100 ka圖5 模擬的不同時(shí)期地下水密度的變化

3.2.2 典型觀測(cè)孔的密度變化

為了對(duì)比不同位置地下水密度隨時(shí)間的變化，在模擬區(qū)域的不同位置、不同深度選取4個(gè)典型觀測(cè)孔，觀測(cè)孔位置見圖4，其中觀測(cè)孔1、2分別位于上層區(qū)域水流系統(tǒng)的補(bǔ)給區(qū)與排泄區(qū)，觀測(cè)孔3、4分別位于下層區(qū)域水流系統(tǒng)的補(bǔ)給區(qū)與排泄區(qū)。不同觀測(cè)孔地下水密度隨時(shí)間的變化見表3。

表3 模擬的觀測(cè)孔地下水密度的變化 kg/m3

位于補(bǔ)給區(qū)的觀測(cè)孔地下水密度保持不變，主要受補(bǔ)給淡水的影響；排泄區(qū)觀測(cè)孔地下水密度持續(xù)增大，且后期的增加速度大于前期的，受蒸發(fā)濃縮作用影響明顯。下層區(qū)域水流系統(tǒng)觀測(cè)孔地下水密度增加幅度大于上層區(qū)域水流系統(tǒng)觀測(cè)孔的，模型運(yùn)行到100 ka時(shí)，位于下層區(qū)域水流系統(tǒng)排泄區(qū)處觀測(cè)孔的地下水密度的增長(zhǎng)值比上層的大4倍多，并且下層區(qū)域水流系統(tǒng)地下水密度的增長(zhǎng)速率也遠(yuǎn)超過(guò)上層區(qū)域水流系統(tǒng)的。

地下水在山前接受補(bǔ)給后，在地形驅(qū)動(dòng)力的作用下，向盆地中心運(yùn)動(dòng)。地下水在補(bǔ)給區(qū)一直保持在淡水狀態(tài)，位于盆地中部的排泄區(qū)受到強(qiáng)烈的蒸發(fā)濃縮作用，溶質(zhì)富集，地下水密度逐漸增大，隨著時(shí)間的推移，地下水由淡水逐漸變?yōu)橄趟?、鹽水，部分地區(qū)最終變?yōu)辂u水。下層區(qū)域水流系統(tǒng)地下水循環(huán)的穿透深度大，地下水徑流路徑長(zhǎng)，排泄區(qū)更接近盆地中心，地下水徑流更加緩慢，有利于鹽分富集和向下部擴(kuò)散，地下水密度以及密度的增長(zhǎng)速率要明顯大于上層的。

3.3 水流驅(qū)動(dòng)力削弱程度的變化

初始時(shí)刻上、下層區(qū)域水流系統(tǒng)H1a與H2a的值如表4所示，用兩者差值表征水流驅(qū)動(dòng)力的大小，則初始時(shí)刻上、下層區(qū)域水流系統(tǒng)地形驅(qū)動(dòng)力的大小分別為69.28、75.58 m。隨著時(shí)間變化，H2a以及δAHC的值均逐漸增大(見表5)，但H1a的水頭基本沒有變化，最終導(dǎo)致H1a到H2a的水流驅(qū)動(dòng)力逐漸減小。

表4 初始條件下區(qū)域水流系統(tǒng)補(bǔ)給區(qū)和排泄區(qū)的平均水頭值 m

受地下水密度差異的影響，上、下層區(qū)域水流系統(tǒng)的δAHC以及E值存在較大差別(見表5)。模型運(yùn)行至10 ka時(shí)，上層區(qū)域水流系統(tǒng)排泄區(qū)的δAHC為2.37 m，E=3.42%。下層區(qū)域水流系統(tǒng)排泄區(qū)的δAHC和E值是上層區(qū)域水流系統(tǒng)的2倍，E=7.12%。運(yùn)行至30～50 ka時(shí)，上層區(qū)域水流系統(tǒng)排泄區(qū)的δAHC和E基本保持不變； 下層區(qū)域水流系統(tǒng)排泄區(qū)的δAHC和E的增長(zhǎng)速率在逐漸增大。模型運(yùn)行至70～100 ka時(shí)，上層區(qū)域水流系統(tǒng)排泄區(qū)的δAHC和E值開始緩慢增大，但E值均小于10%。運(yùn)行至100 ka時(shí)，下層區(qū)域水流系統(tǒng)排泄區(qū)的δAHC為36.17 m，E值高達(dá)47.86%(見圖6)，比上層區(qū)域水流系統(tǒng)的E值高6.5倍，這與高密度流體的分布形態(tài)有很大的相關(guān)性。從圖5中可以看出，下層區(qū)域水流系統(tǒng)的密度增加更為明顯，此時(shí)盆地中部地下水與淡水密度差最大為12%。通常，我國(guó)內(nèi)陸盆地鹽湖分布區(qū)的地下水密度差可達(dá)20%～30%，密度對(duì)地水流動(dòng)力的削弱作用會(huì)更加明顯，區(qū)域水流系統(tǒng)可能受到更大程度的抑制。

表5 區(qū)域水流系統(tǒng)排泄區(qū)平均等價(jià)淡水水頭差變化 m

圖6 上、下層區(qū)域水流系統(tǒng)驅(qū)動(dòng)力削弱程度變化

盆地中部地下水密度增大，對(duì)區(qū)域水流系統(tǒng)水流驅(qū)動(dòng)力造成影響，地下水流線也必將發(fā)生變化。通過(guò)對(duì)比不同時(shí)期地下水流線的徑流距離和深度，可以更直觀地看出密度對(duì)水流驅(qū)動(dòng)力的抑制作用。

圖7所示為從盆地底部出發(fā)的一條流線在不同模擬時(shí)間的位置變化。初始時(shí)刻，流線靠近盆地中部，隨著時(shí)間延長(zhǎng)，逐漸向遠(yuǎn)離盆地中部方向移動(dòng)。模型運(yùn)行至100 ka時(shí)，地下水流線的徑流距離(從補(bǔ)給區(qū)到排泄區(qū)的水平距離)要比初始時(shí)刻縮短了約2 km，徑流深度也減小將近110 m。此時(shí)，區(qū)域水流系統(tǒng)在排泄區(qū)的分布范圍縮小了近1/2，分布范圍被明顯擠壓，區(qū)域水流系統(tǒng)明顯向中間水流系統(tǒng)靠近。

圖7 模擬的不同時(shí)期的流線變化

以上變化表明，盆地中部地下水密度的增大在很大程度上削弱了區(qū)域水流系統(tǒng)的驅(qū)動(dòng)力，改變了地下水徑流路徑，特別是下層區(qū)域水流系統(tǒng)?？梢灶A(yù)見，干旱盆地中部密度差異進(jìn)一步加大的情況下，這種削弱作用也會(huì)進(jìn)一步加強(qiáng)，甚至可能導(dǎo)致區(qū)域水流系統(tǒng)完全消失。

4 結(jié)論

本文中采用平均等價(jià)淡水水頭差分析方法，結(jié)合變密度水流與溶質(zhì)運(yùn)移數(shù)值模擬，研究了干旱盆地中部密度增加對(duì)區(qū)域水流系統(tǒng)驅(qū)動(dòng)力的影響，得到以下主要結(jié)論：

1)理想剖面上的地下水密度變化、單個(gè)觀測(cè)孔的地下水密度變化過(guò)程線分析結(jié)果表明，隨著模擬期的增長(zhǎng)，盆地中部的蒸發(fā)濃縮作用將導(dǎo)致地下水密度增大，等價(jià)淡水水頭差相應(yīng)增大，進(jìn)而導(dǎo)致水流驅(qū)動(dòng)力減小，是影響區(qū)域水流系統(tǒng)的根本原因。

2)盆地中部地下水密度增大，削弱區(qū)域水流系統(tǒng)的水流驅(qū)動(dòng)力，并且對(duì)下層區(qū)域水流系統(tǒng)的削弱程度要遠(yuǎn)高于對(duì)上層區(qū)域水流系統(tǒng)的。模型運(yùn)行至100 ka時(shí)，對(duì)下層區(qū)域水流系統(tǒng)的削弱程度高達(dá)47.86%，比對(duì)上層區(qū)域的削弱程度高約6.5倍。

3)從地下水流線變化上看，地下水密度的增大導(dǎo)致區(qū)域水流系統(tǒng)排泄區(qū)的分布空間被嚴(yán)重?cái)D壓，流線向遠(yuǎn)離盆地中部方向偏移，區(qū)域水流系統(tǒng)的徑流距離和徑流深度都明顯縮小，由此說(shuō)明地下水密度變化對(duì)區(qū)域水流系統(tǒng)影響較大。