基于質(zhì)點運移模型的遼河流域區(qū)域地下水循環(huán)規(guī)律研究

李 昕

(遼寧省遼陽水文局，遼寧 遼陽 111000)

遼河流域是我國東北地區(qū)農(nóng)業(yè)發(fā)展以及糧食生產(chǎn)的重要基地，其地表水資源由于水體污染嚴(yán)重、水資源供需關(guān)系緊張以及水資源浪費等多種因素的綜合作用而日趨嚴(yán)峻，地下水作為該區(qū)域主要供水水源之一對于緩解供水壓力、促進區(qū)域發(fā)展發(fā)揮重要作用。對遼河流域地下水循環(huán)規(guī)律進行研究對于促進區(qū)域水資源可持續(xù)發(fā)展以及科學(xué)規(guī)劃具有重要意義。近年來，國內(nèi)外學(xué)者針對水資源循環(huán)規(guī)律分別采用多個理論和方法從多個角度開展了大量研究，如韓冬梅等[1]采用數(shù)值模擬技術(shù)和地下水流模型對忻州第四系孔隙地下水流平均滯留時間、分布特征進行了分析模擬；孫躍等[2]依據(jù)黃河分布區(qū)域水文地質(zhì)對河南省供水安全問題采用數(shù)值模擬技術(shù)進行研究，并運用局部地下水流預(yù)測模型對10個水源地的地下水開采量進行評價分析；王丹等[3]利用地下水動力學(xué)法對草甸區(qū)剖面地下水進行二維數(shù)值模擬，對地下水運動規(guī)律及影響因素進行揭示和表征；蔣小偉等[4]對地下水模型通過地下水年齡實測法進行矯正，研究表明地下水年齡從補給區(qū)到排放區(qū)整體呈現(xiàn)變老趨勢；地下水年齡在盆地下游垂直方向可產(chǎn)生突變，并以此可對不同級次流動系統(tǒng)進行評判。對遼河流域地下水循環(huán)規(guī)律研究相對較多且主要是采用同位素追蹤法和水化學(xué)法進行探討，而對水動力場的相關(guān)研究主要側(cè)重于利用定性分析，如劉玉等[5]對隴東盆地涇河流域地下水循環(huán)規(guī)律采用水文地球化學(xué)理論和同位素基礎(chǔ)進行研究；徐中華等[6]以隴東盆地白堊系中層地下水為研究對象，采用水文地球化學(xué)理論對其水化學(xué)演化規(guī)律和發(fā)展進行了深入的研究分析；張彥林等[7]對隴東盆地地下水運行狀況采用同位素法和水化學(xué)法進行探討并對其運移規(guī)律以及影響因素進行了探討。

文章采用數(shù)值模擬技術(shù)對遼河流域流速分布特征、水循環(huán)寬度、深度以及地下水運行更新時間進行定量的分析和計算，以期為該流域水資源的科學(xué)規(guī)劃管理以及合理開發(fā)利用提供一定的決策依據(jù)和理論支持。

1 區(qū)域概況

遼河流域位于我國東北地區(qū)西南部，屬于溫帶半濕潤半干旱季風(fēng)氣候，冬季寒冷干燥漫長、夏季炎熱多雨。全年降雨量和徑流量時空分布極不均勻，降雨量由東南向西北方向整體呈降低趨勢，多年平均降雨量為320～860mm，其中62%以上集中在7—8月份，并且多以暴雨或強降雨的形式出現(xiàn)；氣溫在平原地區(qū)較高，山地較低，由南向北整體呈降低趨勢，多年平均氣溫為4～9℃；由南向北方向年蒸發(fā)量依次遞減，并處于982～1650mm范圍。流域內(nèi)地形支離破碎、溝壑縱橫且具有強烈的切割特征，因此為典型的低山丘陵地形，區(qū)域內(nèi)海拔高度由西北向中南部逐漸降低，最高和最低處約為1650m和860m。流域內(nèi)各支流縱橫交錯、蜿蜒曲折，主要支流有老哈河、渾河、太子河、柴河等，研究區(qū)域地質(zhì)結(jié)構(gòu)為白堊系，雖然在流域內(nèi)白堊系分布較為廣泛，但在河谷與溝谷底部或河谷底兩側(cè)多以樹枝狀顯現(xiàn)，受天環(huán)向斜控制，盆地結(jié)構(gòu)多以雁行排列的向斜組成?？砷_采利用的主要含水層位于白堊系含水巖組，而太子河底部含水層巖組和白堊系頂部的砂巖為風(fēng)成砂巖，是主要含水層[8]。渾河底部多以泥巖為主，是該流域的隔水層或弱透水層。地下水多以微淡水、淡水為主，并且受補給條件、含水層厚度以及區(qū)域分布特征作用，該區(qū)域地下水可開采潛力整體處于較低水平。

2 研究方法及內(nèi)容

對遼河流域剖面地下水利用Visual Modflow軟件進行數(shù)值模擬，并且利用向前追蹤例子算法和Modpath模塊生成剖面地下水流線分布，進而建立遼河流域地下水動力場。對遼河流域剖面水頭特征利用地下水動力場進行確定，對其地下水流系統(tǒng)采用Toth地下水流系統(tǒng)理論進行劃分并確定各系統(tǒng)的循環(huán)特征和尺度特征，然后對地下水的平均滯留時間利用Modpath模塊進行計算。

2.1 質(zhì)點運移數(shù)學(xué)模型

Modpath是一種可對任意時刻水質(zhì)點的位置和流線分布進行簡便、快捷的計算和確定的模塊，并可對質(zhì)點向前、向后示蹤流線進行準(zhǔn)確的模擬預(yù)測。Madpath的質(zhì)量平衡方程可對地下水二維穩(wěn)定流動的滲流流速和有效孔隙率進行求解，計算公式如下：

(1)

式中，Vx，Vz—分別為流動流速矢量在坐標(biāo)軸方向的分量，LT-1；w—含水層內(nèi)單位體積水源匯流產(chǎn)生的水量，T-1；n—含水層有效孔隙率，%。

2.2 計算流速

采用單元格的形式對地下水流速度進行計算，即每個單元格對應(yīng)一個水流速度，計算公式如下：

(2)

式中，Vx，Vz—分別為水平和垂直方向的水流流速，LT-1；n—含水層有效孔隙率，%；Qx，Qz—分別為通過水平和垂直方向的水量，L3T-1；Δx，Δy，Δz—分別為單元格在坐標(biāo)軸方向的長度，L。

2.3 徑流劃分標(biāo)準(zhǔn)

首先采用Visual Modflow軟件并結(jié)合每個單元格的徑流速度和剖面地下水流數(shù)值模型提出徑流劃分相關(guān)標(biāo)準(zhǔn)，然后對計算結(jié)果利用surfer10進行差分計算，并完成地下水流速矢量特征的輸出，將其采用色階的模式進行表征?？紤]到遼河流域水流速度基本特征和實際狀況，因此將其徑流劃分為弱徑流、中等徑流和強徑流帶[9]，所對應(yīng)的徑流速度區(qū)間分別為小于0.006m/d、0.06～0.12m/d、大于0.12m/d。

2.4 平均滯留時間

在不同含水層路徑上地下水沿流線方向的平均流速的大小可采用地下水的平均滯留時間進行表征，并采用Modpath放置粒子法即通過對各含水層巖組運動路徑上粒子所消耗的時間進行確定，進而可對地下水的平均滯留時間進行更加直觀、科學(xué)的反映。

3 數(shù)值模型

3.1 剖面概況及概念模型

剖面貫穿遼河流域的東西方向并且在剖面上共布設(shè)18處鉆孔，其中用于環(huán)河含水巖組和含水層研究分析的鉆孔分別為14個和4個。區(qū)域等水位線分布方向整體上與剖面走向呈垂直分布特征，并且白堊系含水巖組厚度在垂直方向上存在顯著差異，并且含水巖組整體上在西部較厚而東部較薄。剖面西部區(qū)域為降雨入滲主要補給區(qū)域，而剖面東西兩側(cè)為主要徑流補給區(qū)域，河道或河道交匯處為水流系統(tǒng)的排泄基準(zhǔn)面并且是地下水排泄的主要排泄區(qū)域[10]。

羅漢洞含水巖組頂板為遼河流域剖面的上邊界，在渾河與太子河交匯處尖滅，并且在上邊界處出現(xiàn)降雨入滲補給，因此該處為流量邊界。側(cè)向徑流補給處往往為剖面兩側(cè)，因此可稱為定流量邊界條件，洛河組含水巖組為流域剖面底部，因此，可稱為隔水邊界。側(cè)向徑流補給和大氣降雨入滲補給為剖面補給的主要來源，地下水向河流的排泄為排泄相的主要來源，采用定水頭的方式對遼河與其他支流交匯處以及沿河道進行處理。含水層的分布特征采用遼河流域剖面水文地質(zhì)進行表征，依據(jù)相關(guān)文獻和研究成果對水文地質(zhì)剖面含水層組的巖性可概化為3層，即洛河、環(huán)河以及羅漢洞含水巖組，可對羅漢洞組、環(huán)河以及洛河含水巖組承壓水分布區(qū)概化為3個、4個和6個滲透系數(shù)分區(qū)，各滲透系數(shù)K1～K13取值為0.005～1.2。通過構(gòu)建數(shù)學(xué)模型進而使得水文地質(zhì)趨于客觀系統(tǒng)的為概念模型的近似處理本質(zhì)，對自然條件下剖面上地下水循環(huán)基本特征進行客觀、準(zhǔn)確的解釋為建立模型的首要目標(biāo)，所以，研究采用非均質(zhì)各向異性穩(wěn)定流對研究區(qū)域模型進行概化，以此提高模型與實際狀況的貼近程度[11]。

3.2 地下水流數(shù)學(xué)模型及求解

可采用下述微分方程對水文地質(zhì)概念模型的定量問題進行求解：

(3)

式中，D，h—分別為剖面滲流區(qū)域和含水層水位標(biāo)高；Kx，Kz—分別為水平方向和垂直方向的滲透系數(shù)，m/d；h1，q(x，z，h)—分別為水頭邊界上水頭分布和流量邊界單寬流量；Γ1，Γ2—分別為水頭邊界和流量邊界。

利用Visula Modflow軟件對模型進行求解，其基本過程為在時間和空間上將研究區(qū)域進行離散化處理，并對每個網(wǎng)格點構(gòu)建均衡方程并聯(lián)立形成一個大型的線性方程組，然后對每個網(wǎng)格的水頭值進行迭代求解。文章將遼河流域剖面劃分為25層、180列和1行的等面積單元格并以此更加準(zhǔn)確、客觀地揭示該區(qū)域的地下水流特征。采用與條件共軛梯度法對水流數(shù)學(xué)模型進行求解，其中內(nèi)部和外部迭代次數(shù)分別為50次，水位變化和殘差收斂標(biāo)準(zhǔn)均為0.01m。

依據(jù)水化學(xué)場特征對模型的水動力場運行結(jié)果進行驗證分析，其中剖面上的水化學(xué)場采用剖面附近12處水井并進行水化學(xué)數(shù)據(jù)分析，并按取樣點位置將所取水樣的水化學(xué)類型進行放置。結(jié)果顯示，地下水陰離子在太子河?xùn)|部區(qū)域類型多變，而且陽離子在排泄區(qū)與過渡區(qū)變化不大，并以鈉離子為主，水化學(xué)特征在整個研究區(qū)域與地下水運移規(guī)律基本保持一致。

4 結(jié)果與分析

4.1 水頭分布特征

剖面模型水頭分布計算結(jié)果顯示遼河流域南部平泉縣附近的水頭為1200～1280m，水頭整體上由南向北方向逐漸降低。該區(qū)域的水頭分布特征為水流由南向北如平泉縣、寧城縣、開魯縣和通遼市等逐漸呈降低趨勢，該區(qū)域?qū)儆谶|河的西部發(fā)源地，共有兩處水源即老哈河和西拉木倫河水源。

4.2 系統(tǒng)劃分

結(jié)合相關(guān)文獻資料將遼河流域剖面依據(jù)地下水流的源、匯關(guān)系以及實際狀況可劃分為區(qū)域循環(huán)、中間循環(huán)和局部循環(huán)三級循環(huán)系統(tǒng)，其中整個剖面內(nèi)的地下水為區(qū)域循環(huán)系統(tǒng)，而在河道、湖泊附近多以中間和局部循環(huán)系統(tǒng)存在。

4.3 尺度特征

由于河流支流眾多，在遼河流域剖面上河流與地下水存在水流交換并形成了多個地下水局部循環(huán)系統(tǒng)，此類循環(huán)系統(tǒng)延伸區(qū)間通常位于5～40km，其發(fā)育深度為20～180m。在老哈河、太子河、西遼河等區(qū)域的局部地下水循環(huán)系統(tǒng)其寬度通常低于10km，其發(fā)育深度約為500～800m；其他主要支流局部循環(huán)系統(tǒng)寬度為10～80km，發(fā)育深度相對較淺，為60～500m。遼河流域北部的碧流河、少冷河、查干木倫河區(qū)間發(fā)育循環(huán)系統(tǒng)，該區(qū)域循環(huán)系統(tǒng)寬度最大為100km，其發(fā)育深度可達到洛河巖組底部最大深度約1200m，各循環(huán)系統(tǒng)的發(fā)育寬度范圍見表1。

表1 各級水流系統(tǒng)發(fā)育范圍

4.4 水循環(huán)特征

在遼河流域剖面上，側(cè)向徑流補給為承壓水的主要補給來源，而大氣降雨入滲補給為流域南部河段的補給源，而北沙河、蘭河、細河作為地下水排泄的主要區(qū)域。分別選取2個過水?dāng)嗝孢M行研究分析，結(jié)果表明地下水流線基本與過水?dāng)嗝嫦啻怪?，而斷面水平徑流在垂直方向流動不顯著。所以，可按照剖面循環(huán)系統(tǒng)對兩個過水?dāng)嗝娼?jīng)過的模型層進行分析，對各單元格的水平流量進行計算和求解，結(jié)果見表2。

表2 遼河剖面三級循環(huán)系統(tǒng)流量分析結(jié)果

研究表明，區(qū)域和中間循環(huán)系統(tǒng)為地下水循環(huán)的主要集中區(qū)域，并且在Z2過水?dāng)嗝嫔暇植亢椭虚g循環(huán)系統(tǒng)的流量占比處于較低水平，而絕大多數(shù)系統(tǒng)流量為區(qū)域循環(huán)系統(tǒng)，其原因主要為洛河組巖層可對Z2斷面層的東側(cè)進行補給，而砂巖的環(huán)河組巖層的補給量相對較低，基本可忽略；結(jié)合區(qū)域?qū)嶋H狀況，環(huán)河組區(qū)域主要為中間循環(huán)系統(tǒng)，而洛河組巖層位于區(qū)域循環(huán)系統(tǒng)。因此，造成區(qū)域循環(huán)系統(tǒng)流量相對中間和局部循環(huán)系統(tǒng)處于較高水平。

4.5 水流分布特征

遼河剖面地下水流平均速度為0.0021m/d，其變化區(qū)間為0.00006～0.012m/d。研究表明：遼河流域的西部區(qū)域?qū)儆谌鯊搅鲙В渌魉俣认鄬^低；渾河區(qū)域地下水流速有一定程度的提升并在洛河組巖層區(qū)達到最大，該區(qū)域?qū)儆趶姀搅骱椭械葟搅鹘唤缣帲魉僭跂|遼河與西遼河區(qū)域逐漸降低為弱徑流帶。整體上看，上部環(huán)河組和下部洛河組分別位于弱徑流帶和中等徑流帶。

4.6 平均滯留時間

研究采用放置粒子法并結(jié)合流速分布特征一級含水層巖組分布狀況對水流運移路徑進行模擬，其中地下水貫穿了3個不同的含水巖層。研究表明，強徑流帶、中等徑流帶和弱徑流帶的循環(huán)更新時間分別為小于4218a，4218～25162a，大于25162a。

5 結(jié)論

剖面模型水頭分布計算結(jié)果顯示，遼河流域南部平泉縣附近的水頭為1200～1280m，水頭整體上由南向北方向逐漸降低。按照區(qū)域特征可將遼河流域劃分為區(qū)域循環(huán)、中間循環(huán)和局部循環(huán)三級循環(huán)系統(tǒng)，其中整個剖面內(nèi)的地下水為區(qū)域循環(huán)系統(tǒng)，而在河道、湖泊附近多以中間和局部循環(huán)系統(tǒng)存在。

在老哈河、太子河、西遼河等區(qū)域的局部地下水循環(huán)系統(tǒng)其寬度通常低于10km，其發(fā)育深度約為500～800m；其他主要支流局部循環(huán)系統(tǒng)寬度為10～80km，發(fā)育深度相對較淺，為60～500m。區(qū)域和中間循環(huán)系統(tǒng)為地下水循環(huán)的主要集中區(qū)域，并且在Z2過水?dāng)嗝嫔暇植亢椭虚g循環(huán)系統(tǒng)的流量占比處于較低水平，而絕大多數(shù)系統(tǒng)流量為區(qū)域循環(huán)系統(tǒng)，其原因主要為洛河組巖層可對Z2斷面層的東側(cè)進行補給，而砂巖的環(huán)河組巖層的補給量相對較低，基本可忽略。

遼河流域北部的碧流河、少冷河、查干木倫河區(qū)間發(fā)育循環(huán)系統(tǒng)，該區(qū)域循環(huán)系統(tǒng)寬度最大為100km，其發(fā)育深度可達到洛河巖組底部，最大深度約1200m；遼河剖面地下水流平均速度為0.0021m/d，其變化區(qū)間為0.00006～0.012m/d；強徑流帶、中等徑流帶和弱徑流帶的循環(huán)更新時間分別為小于4218a，4218～25162a，大于25162a。為提高地下水資源的開發(fā)利用效率和科學(xué)性，可對中等徑流帶的地下水進行合理開采利用，而對弱徑流帶禁止開采或少量開采。